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This project was co-led by Dr. Sheila Mcllraith and Prof. Richard Fikes. Substantial research 
results and published papers describing those results were produced in multiple technology areas, 
including the following: 

Monitoring a Complex Physical System using a Hybrid Dynamic Bayes Net 

The Reverse Water Gas Shift system (RWGS) is a complex physical system designed to produce 
oxygen from the carbon dioxide atmosphere on Mars. If sent to Mars, it would operate without 
human supervision, thus requiring a reliable automated system for monitoring and control. The 
RWGS presents many challenges typical of real-world systems, including: noisy and biased 
sensors, nonlinear behavior, effects that are manifested over different time granularities, and 
unobservability of many important quantities. In this portion of the project, we modeled the 
RWGS using a hybrid (discrete/continuous) Dynamic Bayesian Network (DBN), where the state 
at each time slice contains 33 discrete and 184 continuous variables. We showed how the system 
state can be tracked using probabilistic inference over the model. We investigated how to deal 
with the various challenges presented by the RWGS, and produced a suite of techniques that are 
likely to be useful in a wide range of applications. In particular, we produced a general 
framework for dealing with nonlinear behavior using numerical integration techniques, extending 
the successful Unscented Filter. We also showed how to use a fixed-point computation to deal 
with effects that develop at different time scales, specifically rapid changes occurring during 
slowly changing processes. We tested our model using real data collected from the RWGS, 
demonstrating the feasibility of hybrid DBNs for monitoring complex real-world physical 
systems. 

A Formal Theory of Testing for Dynamical Systems 

Just as actions can have indirect effects on the state of the world, so too can sensing actions have 
indirect effects on an agent’s state of knowledge. In this portion of the project, we investigated 
“what sensing actions tell us”, i.e., what an agent comes to know indirectly from the outcome of a 
sensing action, given knowledge of its actions and state constraints that hold in the world. To this 
end, we developed a formalization of the notion of testing within a dialect of the situation 
calculus that includes knowledge and sensing actions. Realizing this formalization required 
addressing the ramification problem for sensing actions. We formalized simple tests as sensing 
actions. Complex tests are expressed in the logic programming language Golog. We examined 



what it means to perform a test, and how the outcome of a test affects an agent’s state of 
knowledge. Finally, we developed automated reasoning techniques for test generation and 
complex-test verification, and precisely specified restrictions on when the techniques can be used. 
This work is relevant to a number of application domains including diagnostic problem solving, 
natural language understanding, plan recognition, and active vision. 

Diagnosing Hybrid Systems Using a Bayesian Model Selection Approach 

In this portion of the project, we examined the problem of monitoring and diagnosing noisy 
complex dynamical systems that are modeled as hybrid systems; i.e., as systems having 
continuous behavior interleaved by discrete transitions. In particular, we examined continuous 
systems with embedded supervisory controllers that experience abrupt, partial, or full failure of 
component devices. We developed a mathematical formulation of the hybrid monitoring and 
diagnosis task as a Bayesian model tracking and selection problem, and developed a suitable 
tracking algorithm. The nonlinear dynamics of many hybrid systems present challenges to 
probabilistic tracking. Further, probabilistic tracking of a system for the purposes of diagnosis is 
problematic because the models of the system corresponding to failure modes are numerous and 
generally very unlikely. To focus tracking on these unlikely models and to reduce the number of 
potential models under consideration, we exploited logic-based techniques for qualitative model- 
based diagnosis to conjecture a limited initial set of consistent candidate models. We considered 
alternative tracking techniques that are relevant to different classes of hybrid systems, and 
focused specifically on a method for tracking multiple models of nonlinear behavior 
simultaneously using factored sampling and conditional density propagation. A motivating case 
study for this work was the problem of monitoring and diagnosing NASA’s Sprint AERCam, a 
small spherical robotic camera unit with 12 thrusters that enable both linear and rotational 
motion. 

Publications 

The publications listed below describe work that was performed solely or in part in this project. 
The first four of the listed papers are included in this final report. 

U. Lemer, B. Moses, M. Scott, S. Mcllraith, and D. Roller; “Monitoring a Complex Physical 
System using a Hybrid Dynamic Baves Net ”: Proceedings of die Eighteenth Annual 
Conference on Uncertainty in Artificial Intelligence (UA1-2002); Edmonton, Canada; August, 
2002. 

S. Mcllraith and R. Scherl; “What Sensing Tells Us: Towards a Formal Theory of Testing for 
Dynamical Systems ”: Proceedings of the Seventeenth National Conference on Artificial 
Intelligence (AAAI'2000); pp. 483-490; August, 2000. 

S. Mcllraith; “Diagnosing Hybrid Systems: A Bayesian Model Selection Approach ”: Proceedings 
of the Eleventh International Workshop on Principles of Diagnosis (DX'00); Morelia, 

Mexico; June 2000; pp. 140-146. 

S. Mcllraith, G. Biswas, D. Clancy, V. Gupta; “ Hybrid Systems Diagnosis ”: Proceedings of 
Hybrid Systems: Computation and Control; Lecture Notes in Computer Science; Springer- 
Verlag; pgs 282-295; 2000. 

Todd Neller; “Simulation-Based Search for Hybrid System Control and Analysis”; Ph.D. 
Dissertation; Stanford University, Stanford, California, USA; June 2000. 



Monitoring a Complex Physical System using a Hybrid Dynamic Bayes Net 


Uri Lerner 

Computer Science Dept. 
Stanford University 
uri@cs.stanford. edu 


Brooks Moses Mariria Scott Sheila McDraith Daphne Koller 

Mechanical Engr. Dept. Computer Science Dept. Computer Science Dept. Computer Science Dept. 

Stanford University Stanford University Stanford University Stanford University 

bmoses@stanford.edu maricia@cs.stanford.edu sam@ksl.stanford.edu koller@cs.stanford.edu 


Abstract 

The Reverse Water Gas Shift system (RWGS) is a 
complex, physical system designed to produce oxy- 
gen from the carbon dioxide atmosphere on Mars. If 
sent to Mars, it would operate without human super- 
vision, thus requiring a reliable automated system for 
monitoring and control. The RWGS presents many 
challenges typical of real-world systems, including: 
noisy and biased sensors, nonlinear behavior, effects 
that are manifested over different time granularities, 
and unobservability of many important quantities. In 
this paper we model the RWGS using a hybrid (dis- 
crete/c ontinuous) Dynamic Bayesian Network (DBN), 
where the state at each time slice contains 33 discrete 
and 184 continuous variables. We show how the sys- 
tem state can be tracked using probabilistic inference 
over the model We discuss how to deal with the var- 
ious challenges presented by die RWGS, providing a 
suite of techniques that are likely to be useful in a 
wide range of applications. In particular, we describe 
a general framework for dealing with nonlinear behav- 
ior using numerical integration techniques, extending 
the successful Unscented Filter. We also show how 
to use a fixed-point computation to deal with effects 
that develop at different time scales, specifically rapid 
changes occurring during slowly changing processes. 

We test our model using real data collected from the 
RWGS, demonstrating the feasibility of hybrid DBNs 
for monitoring complex real-world physical systems. 

1 Introduction 

The Reverse Water Gas Shift System (RWGS) shown in 
Fig. 1 is a complex physical system designed and con- 
structed at NASA’s Kennedy Space Center to produce oxy- 
gen from carbon dioxide. NASA foresees a number of pos- 
sible uses for the RWGS, including producing oxygen from 
the atmosphere on Mars and converting carbon dioxide to 
oxygen within closed human living quarters. 

In a manned Mars mission, the RWGS would operate 
for 500 or more days without human intervention [Larson 
and Goodrich, 2000]. This level of autonomy requires the 
development of robust and adaptive software for fault diag- 
nosis and control. In this paper, we focus on two key sub- 
tasks — monitoring and prediction. Monitoring, or track- 
ing the current state of the system, is a crucial component 



Figure 1: The Prototype RWGS System 


of the control system. Prediction of the system’s expected 
behavior is a basic tool in fault diagnosis — discrepancies 
between the predicted and the actual behavior of the system 
may indicate the presence of faults. 

The RWGS presents a number of significant modeling 
and algorithmic challenges. From a modeling perspec- 
tive, the system is very complex, and contains many sub- 
tle phenomena that are difficult to model accurately. Var- 
ious phenomena in the system manifest themselves over 
dramatically different time scales, ranging from pressure 
waves that propagate on a time scale of milliseconds to 
slow changes such as gas composition that take hours to 
evolve. From a tracking perspective, the system dynamics 
are complex and highly nonlinear. Furthermore, the sen- 
sors give only a limited view of the system state. Some key 
quantities of the system are not measured, and the available 
sensors are noisy and biased, with both the noise level and 
the bias varying with the system state. 

In this paper we model the RWGS using a hybrid (dis- 
crete/continuous) Dynamic Bayesian Network (DBN), and 
show how the system state can be tracked using probabilis- 
tic inference over the model. We focus on the continuous 
part of the model, assuming all the discrete variables are 
known. We discuss how to deal with the various challenges 
presented by the RWGS, both in terms of modeling and in 
terms of inference. We provide a suite of techniques that 
are likely to be useful in a wide range of applications, in- 


eluding the ease where the discrete variables are not ob- 
served. 

Perhaps the most interesting modeling problem pre- 
sented by the RWGS is the issue of different time granu- 
larities. A naive solution is to discretize time at the finest 
granularity. Unfortunately, this approach is generally in- 
feasible both because of the computational burden and be- 
cause the number of observations is effectively reduced to 
one for every few thousand time steps, leading to serious 
inaccuracies. Instead, we take the approach of modeling 
the system at the time granularity of the observations. We 
show how to deal with the almost instantaneous changes 
relative to our time discretization by modeling a part of our 
system as a set of fixed-point equations. 

For the inference task, we provide some new insights 
into the problem of tracking nonlinear systems. This task 
is commonly performed using the Extended Kalman Fil- 
ter (EKF) [Bar-Shalom et al L, 2001] or the simpler and 
more accurate Unscented Filter (UF) [Julier and Uhhnann, 
1997]. We view the problem as a numerical integration 
problem and demonstrate that the UF is an instance of a 
numerical integration technique. More importantly, our ap- 
proach naturally leads to important generalizations of the 
UF: We show how to take advantage of the structure of the 
DBN and present a spectrum of filters, trading off accuracy 
with computational effort. 

We tested our model using real data collected from the 
RWGS prototype system. Our results demonstrate the po- 
tential of using hybrid DBNs as a monitoring tool for com- 
plex real-world physical systems. 

2 Preliminaries 

In this paper, we characterize physical systems as discrete- 
time stochastic processes. System behavior is described in 
terms of a system state which evolves stochastically at dis- 
crete time steps / = 0-1.2 We assume that the system 

is Markovian and stationary , i.e., the state of the system 
at time / -f 1 only depends on its state at time /, and the 
probabilistic dependencies are the same for all / . 

The system state is modeled by a set of random vari- 
ables .V - (.Vi X v ] . We partition the state vari- 

ables .V into a set of evidence (observed) variables, E , and 
a set of hidden (unobserved) variables, H. Physical sys- 
tems commonly comprise both continuous quantities (e.g., 
flows, pressures, gas compositions) and discrete quanti- 
ties (e.g., valve open/closed, compressor on/off). Conse- 
quently, we model such systems as hybrid systems, with X 
comprising both discrete and continuous variables. 

We model the process dynamics of our system using a 
Dynamic Bayesian Network (DBN) [Dean and Kanazawa, 
1989]. A DBN is represented as a Bayes Net fragment 
called a 2TBN, which defines the transition model P{ X* | 
X) where X f — f.Yi — , V fr , ] denotes the variables at 
time / -f 1 and X = (A'i , . . . , .V,/] denotes some subset of 
the variables at time / which are persistent ; in that their val- 
ues directly influence the next state. More formally, a DBN 


is a directed acyclic graph, whose nodes are random vari- 
ables in two consecutive time slices, X and X' . The edges 
in the graph denote direct probabilistic influence between 
the parents and their child. For every variable V' at time 
/ -f 1 we denote its parents as Par( .Y') C X U X'. Each 
.V' is also annotated with a Conditional Probability Dis- 
tribution (CPD), that defines the local probability model 
P(X l | Par( V' )). In our hybrid model, discrete nodes do 
not have continuous nodes as parents. 

The tracking problem in DBNs is to find the belief state 

distribution Bel(X') P{X f \ r* 1 , r'), where X 1 typ- 
ically consists of the persistent variables X at time /, and 
r 1 . c 1 are the evidence variables from time 1 to time I . 
The belief state summarizes our beliefs about the state of 
the system at time /, given the observations from time 1 
to time / . As such, it makes current and future predictions 
independent of past data. The tracking algorithm is an it- 
erative process that propagates the belief state. We start 
with the belief state at time /, Bel(X') and perform three 
steps. We first compute P(X\ X' +1 | r 1 , ...,r / ) as the 
product Bel(X i )/ , (X /+1 | X'). Next we marginalize out 
X' resulting in a distribution over X /+1 . Finally, we con- 
dition on and the result is the belief state at / -f 1, 
Bel(X /+1 ). 

Linear models are an important class of DBNs. In a 
linear model, all the variables in X are continuous and 
all the dependencies are linear with some added Gaussian 
noise. More precisely, if a node X has parents Yi , .... Y* 

then P{ X | Yi H) = ]T* =1 + V, where the ?/ r ’s 

are constants and V has a normal distribution A'i/t.rr 2 ). 
In a dynamic linear model, tracking can be done using a 
Kalman filter [Kalman, I960], where the belief state is rep- 
resented parametrically as a multivariate Gaussian in terms 
of the mean vector and the covariance matrix. Kalman fil- 
ters therefore allow a compact belief state representation, 
which can be propagated in polynomial time and space. 

When the dependencies in the model are nonlinear, the 
resulting distributions are generally non-Gaussian and can- 
not be represented in closed form. Consequently, the belief 
state is generally approximated as a multivariate Gaussian 
that preserves the first two moments of the true distribu- 
tion. The traditional method for doing this approximation 
is using an Extended Kalman Filter (EKF) [Bar-Shalom et 
al , 2001]. Assume that X* — /(X), where / is some 
nonlinear function and X ^ ,V(//. Y). Note that we can 
always assume that f is deterministic: If the dependency 
between X and X 1 is stochastic we can treat the stochas- 
ticity as extra random variables that f takes as arguments. 
The EKF finds a linear approximation to / around the mean 
of X, i.e., we approximate / using the first-order Taylor 
series expansion around //. The result is the linear func- 
tion /(X) » /(/#) + Yf\ ft (X - /#), where Yf\ ft is the 
gradient of / evaluated at p} 


l A second-order EKF approximation exists, but its increased 
complexity tends to limit its use. 



The EKF has two serious disadvantages. The first is its 
inaccuracy — the EKF is accurate only if the second and 
higher-order terms in the Taylor series expansion are neg- 
ligible. In many practical situations, this is not the case 
and using the EKF leads to a poor approximation. The 
second disadvantage is the need to compute the gradient. 
Some nonlinear functions may not be differentiable (e.g., 
the max function), preventing the use of an EKF. Even 
if the function is differentiable, computing the derivatives 
may be hard if the function is represented as a black box 
rather than in some analytical form. 

The Unscented Filter (UF) [Julier and Uhhnann, 1997] 
provides an alternative approach to tracking nonlinear be- 
havior. As with the EKF, the UF assumes that X r = /( X ) 
and X — A'(/#, £). The UF works by deterministically 
choosing 2d ■+ 1 points *o, * 2./* where arn = /# and the 
other points are symmetric around // (die actual points de- 
pend on ^). Associated with each point is a weight tiy. 
The UF computes x J =■ /( 37 ) for i = 0, 1 , . . . , 2 d 9 result- 
ing in 2d + 1 points in 1R T7? , from which it estimates the first 
two moments of X 1 as a weighted average of the In 

particular, the mean F[X'] is approximated as Y/ito «’«*»{ • 

The UF has several significant advantages over the EKF. 
First, it is easier to implement and use than the EKF — no 
derivatives need be computed, and the function / is simply 
applied to 2d + 1 points. Second, despite its simplicity, the 
UF is more accurate than the EKF: The UF is a third-order 
approximation, i.e., inaccuracies are induced only by terms 
of degree four or more in the Taylor series expansion. Fi- 
nally, instead of just ignoring the higher-order terms, the 
UF can account for some of their effects, by tuning a pa- 
rameter used in the point selection. As shown in [Julier and 
Uhlmann, 1997], the UF can be extremely accurate, even in 
cases where the EKF leads to a poor approximation. 

3 The RWGS System 

The purpose of the RWGS is to decompose carbon diox- 
ide (C0 2 ) (abundant on Mars) into oxygen (O2) and carbon 
monoxide (CO). The system, shown in Fig. 2(a) [Goodrich, 
2002], comprises two loops: a gas loop that converts CO? 
and hydrogen (H 2 ) into H 2 0 and CO, and a water loop that 
electrolyzes the H 2 0 to produce 0 2 and H 2 . Under normal 
operation, C0 2 at line (1) is combined with H 2 returned 
from the electrolyzer via line (12), and a mixture of C0 2 , 
H 2 , and CO from the reactor recycle line (11). This mixture 
enters a catalyzed reactor (3) heated to I00°C. Approxi- 
mately 10% of the COo and H*> react to form CO and H 2 0: 
C0 7 + IF =±CO - h IUO 

The H 2 0 is condensed at (4) and is stored in a tank (5). The 
remaining gas mixture passes through a separation mem- 
brane (9), which sends a fraction of the CO to the vent (13) 
while directing the remaining mixture into the recycle line 
( 1 1). A compressor ( 10) is used to maintain the necessary 
pressure differential across the membrane. In the water 
loop, the H 2 0 in tank (5) has some C0 2 dissolved in it, 
which would be detrimental to the electrolyzation process. 


To remedy this, the H?0 is pumped into a second tank (6), 
and has Ho bubbled through it to purge the COo. From 
there, the H 2 0 is pumped into the electrolyzer (8), which 
separates a portion of it into 0 2 and H 2 . The H 2 re-enters 
the gas loop via ( 12), while the remaining H 2 0, along with 
the 0 2 , goes into tank (7), where the mixture is cooled and 
separated. The HoO returns to the electrolyzer, while the 
0 2 leaves the system through ( 14). 

In addition to its normal operating mode, the system 
may operate without the electrolyzer and water pumps. In 
this mode, the Ho for the reaction is supplied by a supply 
line (15) paralleling the C0 2 supply line. This option is not 
feasible for operation on Mars, but has proven useful for 
testing the physical system while under development. 

The RWGS is an interconnected nonlinear system 
where the various components influence each other in com- 
plicated and sometimes unexpected ways. For example, 
during runs without the electrolyzer, it is necessary to 
empty the water tank (5) periodically, to prevent water from 
accumulating and eventually overflowing the tank. This 
causes the gases in the tank to expand, and thus creates a 
significant and sudden pressure drop, which affects the flow 
throughout the whole system. This phenomenon is demon- 
strated in Fig. 2(b), taken from [Whitlow, 2001]. The graph 
shows the flow through the CO vent ( 13) as it evolves over 
time — the spikes correspond to emptying the water tank. 

A challenging property of the RWGS is that phenomena 
in the system manifest themselves over at least three differ- 
ent time scales. Pressure waves in the RWGS propagate 
essentially instantaneously (at the speed of sound). Gases 
flow around the gas loop on the order of seconds. Finally, 
gas compositions in the gas loop take on the order of hours 
to reach a steady state. Meanwhile, the sensors collect data 
at a sampling rate of one second. 

An additional challenge of the RWGS is its sensitivity 
and unidentifiability, i.e., parts of the system state are very 
sensitive to input paramaters and are not directly measured. 
For example, the H 2 and C0 2 compositions in the gas loop 
cannot be practically measured. However, the balance be- 
tween these compositions is almost neutrally stable, thus 
a small shift in the input conditions or the membrane be- 
havior will cause the balance to gradually drift to a signifi- 
cantly different value. 

As in any real system, the RWGS sensors do not record 
the underlying state exactly. In addition to some impor- 
tant quantities, such as the gas compositions, which are not 
measured at all, the existing sensors are noisy and biased. 
The noise level of the sensors depends on many factors and 
can change over time. An example is shown in Fig. 2(c), 
where the difference in the readings of the pressure sensors 
P3 and Pa (both located at (2) in Fig. 2(a)) is plotted over 
time. The main reason for the noise in time steps 0-800 
is the physical proximity of the sensors to the compressor 
that sends pressure waves throughout the system. Since 
the sensors are not synchronized with the compressor, they 
take measurements at various phases of the pressure waves 
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Figure 2: (a) The RWGS Schematic, (b) Effects of emptying a water tank, (c) Pressure difference between P* and P 4 . 


and thus measure significantly different values. After 796 
seconds the compressor shuts down and the noise level de- 
creases dramatically. 2 More interestingly, we note that 
the two sensors are placed very close together and thus the 
average difference should be zero. However, as the plot 
demonstrates, this is not the case, indicating that the sen- 
sors are not well calibrated and some bias is present. Fur- 
thermore this bias depends on the system state, as shown by 
the change in the average difference when the compressor 
shuts off. 

4 Modeling the RWGS 

We model the RWGS using a hybrid DBN, as described in 
Section 2. The 2TBN has 293 nodes, 227 of which are con- 
tinuous. Currently the discrete variables in the model are 
all known and correspond to computer-controlled switches 
and sensor faults. The continuous variables in our model 
capture the continuous-valued elements of our system (e.g., 
pressure at various points in the system, flow rates, tem- 
peratures, gas composition, etc.). Of the 227 continuous 
nodes, 43 represent the time / belief state X and 1 84 repre- 
sent the variables X' at time i + 1 . Of the latter, 43 variables 
are belief state variables for / -f 1 , 72 variables are encap- 
sulated variables, as discussed in Section 5.4, and the rest 
are either sensor variables or transient variables. 

When constructing the model, we used four techniques 
for parameter estimation. Some of the parameters were 
known physical constants or system properties. Of the em- 

2 The sensor’s noise is literally noise that can be heard — the 
pressure waves are the sound waves generated by the compressor. 


pin cal parameters, many came from physical models. The 
others (specifically, some parameters for the compressor, 
the separation membrane and the overall system pressure 
changes) were determined using common equations that 
model the particular system behavior. All the variables in 
these equations were directly observed in the data, and thus 
we could use least-squares techniques to find the best fit for 
the parameters. The remaining parameters were estimated 
using prior knowledge of the domain. 

4.1 Sensor Modeling 

As discussed in Section 3, one of the challenges we address 
in modeling the RWGS is dealing with noisy and biased 
sensors. We deal with noisy sensors in the obvious way: 
by increasing the variance of the predicted measurement 
values to match the noise level in the data. 

Sensor biases present a more interesting modeling prob- 
lem. The biases are not easily modeled using a simple pa- 
rameter since they are unknown and can drift over time. In- 
stead, we address the problem by adding hidden variables 
to the belief state that model the different biases of the sen- 
sors. Biases start with zero mean and a reasonably large 
variance and persist over time, i.e., Bias' =r Bias' + V, 
where \ ' represents white noise with a relatively small vari- 
ance, allowing for some amount of drift to occur over time. 

This idea works quite well, but it tends to overfit the 
data: By letting the bias account for every discrepancy be- 
tween the model predictions and the actual sensor measure- 
ments, the tracking algorithm might settle in an incorrect 
steady state. To fix the problem we must make sure that 
the model biases reflect true sensor biases — biases should 



be kept as small as possible and allowed to grow only if 
there is a real reason for that. We implement this idea by 
introducing a contraction factor -> < 1 (empirically set to 
be 0.07) into the bias formula: Bias'* 1 = -) * Bias' *f V\ 
Thus, biases tend to go to zero unless doing so introduces 
a systematic discrepancy with the predicted system state. 

4.2 Sensitivity and Unidentifiability 

Recall that the equations governing the H 2 /CO 2 balance in 
the gas loop are sensitive to slight variations in the physical 
parameters. Thus even using the most exact form of these 
equations in the model will result in (at least) the same level 
of sensitivity — both to variations in the physical parame- 
ters, and inherent errors in the parameters. Moreover, the 
model value is also sensitive to model effects such as cal- 
culation errors and sensor errors that do not affect the real 
value. We therefore use equations for the IWCO? balance 
that contain an intentionally non-physical component — a 
stabilizing teim — that reduces the sensitivity. This term 
drives the balance to a pre-detennined point, which in this 
case is our expected value for the balance. The magnitude 
of this term is manually adjusted to provide an optimum 
tradeoff between physical accuracy and model stability. 

43 Differing Time Scales 

As described in Section 3, we must deal with differing time 
scales in modeling the RWGS. The naive solution to this 
problem is to model the DBN at a very fine time granularity. 
However, it is completely impractical to model the behav- 
ior of the pressure waves using a discretized-time model. 
To do so would require time steps three orders of magni- 
tude smaller than the time between measurements, which 
is a significant waste of resources. Furthermore, it would 
require a much more complete description of the system 
than is practical, and tracking the slowly-evolving aspects 
of the system with a step size many orders of magnitude be- 
low their time scale would allow substantial errors to build 
up. 

Thus, we approximate the pressure waves as occurring 
instantaneously and instead of modeling their transient be- 
havior, we model the quasi-steady-state results at each time 
step after they have reached an equilibrium. The equa- 
tions in this case are substantially simpler, and require far 
fewer empirical constants. The difficulty, however, is that 
these equations must be solved simultaneously; a change 
in any part of the system will affect all of the other parts. 
These equations include both the compressor equation and 
an approximation to the membrane equations developed 
in [Whitlow, 2001]; thus, they are fairly large and nonlin- 
ear, and no direct simultaneous solution form exists. In- 
stead, we use these equations to create a new equation that 
converges to a fixed point solution. 

We must insert this fixed-point equation into a (nonlin- 
ear) CPD to use it in our DBN model of the RWGS. The 
equation solves for the five model variables Z that account 
for the flows and pressure of file gas loop. In order to solve 


for all five variables, their eight parents must also be present 
in the CPD. Hence, we have a vector CPD for Z whose def- 
inition is essentially procedural: given a value of the eight 
parents it executes an iterative fixed-point computation un- 
til convergence, and outputs the values Z. 

5 Tracking In Nonlinear Systems 

In this section, we address the problem of inference, fo- 
cusing on tracking in complex nonlinear systems, such as 
the RWGS. In these models, the probabilistic dependen- 
cies, including sensors, can be either linear or nonlinear 
functions with Gaussian noise. We restrict our attention to 
the task of tracking the continuous state, assuming all the 
discrete values are known. Note that although the results in 
this section are presented in terms of dynamical systems, 
the analysis also applies to probabilistic inference in static 
nonlinear Bayes nets. 

5.1 Exploiting DBN Structure 

Recall the setup from Section 2: We have a Gaussian belief 
state Bel(X) where X £ 1R / and a 2TBN representing 
P(X ' | X) as a deterministic function X* — f(X). Our 
goal is to find an approximation of P{ X f ) as a multivariate 
Gaussian. The classical approach, used in the EKF and 
the UF, is to find the entire distribution P{ X') directly by 
treating / as a function from 1R' 7 to JR ,r \ An alternative 
approach is to decompose / by defining AT — //(TV) for 
7 = 1 , . . . . in f where Y, ~ Par( AT). In most practical 
cases the /, ’s have a lower dimension than /; as we shall 
see, this reduction in the dimension lets us approximate the 
resulting distribution more accurately and efficiently. 

As discussed in Section 2, the first step in the be- 
lief state propagation process is to compute a multivari- 
ate Gaussian over (X. X*}. We begin with our Gaus- 
sian Bel(X), and add the variables from X' one at a 
time, using the procedure described in Section 5.2. The 
key insight is that, as X- is conditionally independent of 
fX — Yi, -V] AT., ] given Y ti it suffices to approx- 

imate the Gaussian P(TV AT). We can then compute 

rix. av..,at)= HX,.v; at_,)p(at | y.), 

which, for Gaussians, can be accomplished using simple 
linear algebra operations. 

A more difficult case arises when the DBN contains not 
only inter-temporal edges from X to X f , but also intra- 
temporal edges between X f variables. In this case we 
sort the variables Xj in topological order, and gradually 

build up the joint distribution P(X , .V] AT). The 

topological order ensures that when we need to compute 
P(TV AT), we have already computed a Gaussian over 

Yi C XUf .VJ AT., } . This approach, however, may 

introduce some new inaccuracies, because we now also use 
a Gaussian approximation for the distribution of the rele- 
vant variables from {.Y{ , AT_i ] . 

Even in cases where we introduce extra inaccuracies, 
this method is often superior to the UF. The reason is that, 
by reducing the dimension of the functions involved, we 



can use more accurate techniques to approximate the first 
two moments of the variables in X' with the same compu- 
tational resources. In general, there is a tradeoff between 
the superior precision we achieve for each variable and the 
potential for extra inaccuracies we introduce. The extra in- 
accuracies depend on the quality of our Gaussian approxi- 
mation for P(X. X\ .Y^), and on the extent of the 

nonlinearity of the dependencies within X If the depen- 
dence of Y' on f X \ , . . . , Y ■ _ 1 ] is linear, then there are 
no extra errors introduced: In this case the first two mo- 
ments of X} are only influenced by the first two moments 
of f Y' , . . . , .Y;_! } which can be captured correctly by our 
Gaussian approximation. It is somewhat reassuring that the 
better our approximation of P(X f ) as a Gaussian is, the 
less significant the extra errors we introduce are, as the en- 
tire framework is based on the assumption that P( X 1 ) can 
be well approximated by a Gaussian. 

5.2 Numerical Integration 

We now turn our attention to the task of approximating 
P{Y ,Y f 0 as a multivariate Gaussian. To simplify our 
notation, let X be a variable which is a nonlinear func- 
tion of its parents Y — Vi,... .V?, i.e., X - f{Y\ 
but the ensuing discussion also holds for the vectoT case of 
X - f{Y). We assume that P(Y) is a known multivariate 
Gaussian, and the goal is to find a Gaussian approximation 
for P(Y, .Y). It suffices to compute the first two moments: 


m = 

j P{Y)J(Y)dY 

( 1 ) 

nx*} = 

f P{Y)f-{Y)dY 

( 2 ) 

n[\Yj] = 

j P{Y)f{Y)YjdY 

( 3 ) 


Note that the integrals only involve the direct parents 
of Y, significantly reducing their dimension. We can ef- 
fectively compute these integrals using a version of the 
Gaussian Quadrature method called the Exact Monomial 
rules [Davis and Rabinovitz, 1984]. Generally speaking, 
Gaussian Quadrature approximates integrals using a for- 
mula of the form: 

f v 

I W(x)f(x)dx fis ^ «•;/(*;) 

J - 1 

where IT 7 ( x) is a known function (a Gaussian in our case). 
The points x ; and weights vj are carefully chosen to en- 
sure that this approximation is exact for any polynomial / 
whose degree is at most p. The degree p is called the pre- 
cision of the approximation. 

Finding a set of points with a minimal size V for some 
precision p is not a trivial task. In the simple form of Gaus- 
sian Quadrature, we choose points in one dimension and 
use them to create a grid of points in IR J with the obvious 
disadvantage that V grows exponentially with d. Fortu- 
nately, we can do better. In [McNamee and Stenger, 1967] 


0.30 

Optimal 

or 


. Intagmfcn 

^ From sampte 


PrtiCttonS 

32* 

/• •* . 


PracMon S ¥ Y\ 

.. * s 

0.20 

, 1" \ 


* ♦ . 


/ / \ \ 

326 

, ,7i 


/ / v - \ 



0.10 

/ jl \\ 

ZM 



/J/ 


-V' i • 

-4 0 4 

(a) 

SOO 3L0 306 40 

(b) 


Figure 3: (a) Density estimates for Y = y/Yp -Mo 2 . (b) 
Random samples from the RWGS network for the flow at 
point (16) and the pressure at point (2), and Gaussian esti- 
mates for the distribution. 

a general method is presented for .V - O (i^-) and pre- 
cision p - 2k + 1 (d is the dimension of the integral, in 
our case \Y\). In particular, rules are presorted for 2d + 1 
points with precision 3, 2d 2 + 1 points with precision 5 
and + §</+! points with precision 7. The precision 3 
rule is exactly the rule used for the Unscented Filter: It has 
exactly the same 2d + 1 points and weights. 

This view of the Unscented Filter has immediate prac- 
tical consequences: we can trade off between the accuracy 
of the computation and its computational requirements. For 
example, if we are interested in a more precise filter than 
the Unscented Filter and are willing to evaluate the func- 
tion at f>(d 2 ) points then we can use the exact monomial 
rule of precision 5. Depending on the function, this may 
represent a significant gain in accuracy. 

As a simple example we consider the nonlinear function 
.V = yVf + V, 2 where P(V,) = A'(2,-l) and P(V 2 | 
Vi) = A'(0.5Vi - 1,3) (note that both V, andV 2 have die 
same variance I). Fig. 3(a) shows various estimates for the 
probability of X. The optimal estimate is the best Gaussian 
approximation for the distribution of A' computed using a 
very exact numerical integration rule. We can see that the 
exact monomial rules of precisions 3 and 5 provide a much 
better estimate than EKF, where the precision 5 rule leads 
to a more accurate estimate than the precision 3 rule. 

53 Inaccuracies In the Approximation 

Unfortunately, approximating P(Y, X) using numerical 
integration can lead to covariance matrices that are not 
semi-positive definite, and hence illegal. One simple ap- 
proach to this problem is to use a more accurate integration 
rule, although the problem may persist. An alternative is to 
find the “closest” positive definite covariance matrix. We 
cast this problem as a convex optimization problem follow- 
ing [Boyd and Vandenberghe, 2003]. 

Consider once again the problem of approximating 
P[Y , Y) as a multivariate Gaussian, where Y is a non- 
linear function of its parents Y, i.e., X = f(Y\ and 
Y ~ A'(/<v» -rr). Let V denote the estimated covari- 



ance matrix for P(Y, .V): 



If n and f lead to a matrix 31 that is not positive definite, 
then we need to find the closest u and r to fi and r, such 
that Y is positive definite. Given that 3yy is already pos- 
itive definite, 31 is positive definite iff t — » T Ey y M > 0. 
Thus, we can formalize our problem as follows: 

Minimize || n - n || 2 +(r - r)~ (4) 

Subject to n r 3yyU — r + c < 0 (5) 

where c is some small positive number. Since both Eq. 4 
and Eq. 5 are convex we can solve this problem by form- 
ing the Lagrangian and solving the dual problem. We set 
the partial derivatives of it and r to zero and plug the result 
into Eq. 5. We get an equation over the Lagrangian multi- 
plier which can be solved easily as it involves a monotonic 
function. We omit details for lack of space. 

Our analysis treats the elements in n and r directly, 
but in fact these elements are not independent since i/,* — 
n[Y,X] - // r , n[.X] and r = r[X~] - 77[.Y] 3 * * * * * . It is desir- 
able to use this relation in Eq. 4 and Eq. 5 and represent the 
dependency between the various elements (e.g., a change in 
/7[.Y] may fix many of the problems simultaneously). Un- 
fortunately, because of the tenn 77[.Y] 2 the problem is no 
longer convex. Nonetheless, we can approximate the prob- 
lem as convex (e.g., by replacing 77 [.Y] 2 by the best current 
estimate), solve it and iterate. Again, we defer details to an 
extended version of this paper. 

5.4 Encapsulated Variables 

Just as we can use the DBN structure to decompose the de- 
pendency between X 9 and X , in many cases we can further 
decompose the dependency Y — /( Y). For example, as- 
sume that f(Y) — g{g\ (Y i ). g?(Y?)), where Yi , Y? C 
Y ? Instead of directly evaluating the Gaussian over 
(Y, .Yj we can define two extra variables: Ti = g\ ( Yy ) 
and T 2 = < 7 ?(Y 2 ). We first approximate P(Yi.Ti) as 
a Gaussian and use it to find a Gaussian over {Y,T|}. 
Next we approximate P[ Y 2 , To) as a Gaussian and from it 
P( Y. Tt , T 2 ). Finally, we approximate P(T i , To, Y) as a 
Gaussian and use it to find the Gaussian approximation for 
P( Y. Ti , To, .Y). The same accuracy tradeoffs that were 
discussed in the context of X 1 — J(X) apply here: by 
reducing the dimension of the integrals we can solve each 
one more accurately, but may introduce further errors if the 
interaction between the extra variables is nonlinear. 

3 E.g., flow sensors give different results depending on the gas 

type. Assuming we have random variables representing the total 

flow and the compositions of the different gases in it, and 

may each be a product of one of the gas compositions and the 

flow, thus representing the net flow of a certain gas. The func- 

tion g would be a weighted sum of these flows where the weights 

correspond to the sensor’s response for the different gases. 



Figure 4: Comparison with particle filtering on simulated 
data, showing the means and error bars of two standard de- 
viations for our algorithm and particle filtering. The Y axis 
represents time, and the 3 ’ axis the percentage of H 2 in the 
flow at point (16). To increase readability, we shift the es- 
timates generated by our algorithm by 0.1 to the left and 
those generated by particle filtering by 0.1 to the right. 

In principle, one could add T f and T 2 to the DBN and 
treat them as regular variables. However, doing so makes 
these variables part of X\ and thereby increases the al- 
gorithm’s space complexity, which is 0{ [X'l 2 ) (for repre- 
senting the covariance matrix of P[ X 1 )). It is better to treat 
the extra variables as local variables encapsulated within 
the CPD and unknown to the rest of the network. After 
computing the Gaussian approximation for the CPD vari- 
ables, we simply marginalize over the encapsulated ones. 
This approach is similar to the local computations in an 
OOBN model [Koller and Pfeffer, 1997], where some of 
the CPD variables are encapsulated within the CPD. 

6 Experimental Results 

In this section we present results from a set of experiments 
that test the efficacy and robustness of our model and track- 
ing algorithm. Our computational model of the RWGS con- 
tains all of the components needed to monitor the full op- 
eration of the physical system, although data provided to 
date by KSC is for the reduced-operation mode with only 
the gas loop component operational. Our experiments were 
run on a Pentium III 700MHz. 

We tested our algorithm with both real data and simu- 
lated data that was generated from our model. Although 
running with real data is the real test for our approach, run- 
ning with simulated data is also of interest. The reason is 




that there are two sources of errors when using real data: 
model inaccuracies and errors induced by the algorithm. 
When using simulated data, only errors of the second type 
are present and we can better test the performance of the 
algorithm. 

6.1 Results on Simulated Data 

We first tested whether the belief state could be well ap- 
proximated as a Gaussian and whether our particular ap- 
proximation was a good one. To do so, we generated a set 
of samples from the model. We did not introduce any evi- 
dence so the samples were indeed sampled from the correct 
joint distribution. In Fig. 3(b) we show the results for two 
particular variables: the flow at point ( 1 6) and the pressure 
at point (2) (these variables were chosen because of their 
dependency on the non-linear CPD of the membrane; other 
variables produced similar results). The samples appear to 
be drawn from a distribution that is either a Gaussian or 
close to one. Furthermore, our estimate for the joint distri- 
bution (depicted by the contours for one and two standard 
deviations) is very close to the Gaussian that was estimated 
directly from the samples. Thus, it is reasonable to expect 
that our techniques will lead to good approximations of the 
belief state. 

Next, we generated a trajectory of 500 time steps from 
our model and tested our algorithm on it. We compared 
our results with the particle filtering algorithm [Gordon et 
al . , 1993], which approximates the belief state as a set of 
weighted samples where the weights of the samples corre- 
spond to the likelihood of die evidence given the sample. 
Our algorithm took 20ms per time step, which included 
computing the Gaussian approximation to the belief state, 
with numerical integration when necessary, and condition- 
ing on the evidence. In comparison, generating a sample 
using particle filtering took 1 .5ms. Thus, one step of our al- 
gorithm took as much time as generating 13 samples. How- 
ever, with just 13 samples particle filtering performed ex- 
tremely poorly and therefore in our experiments we used 
10,000 samples at every time step, giving particle filtering 
a somewhat unfair advantage. 

Fig. 4 shows the estimates for the percentage of Ho in 
the flow at point ( 16) that were computed by our algorithm 
and by particle filtering, as well as the actual value (known 
from the simulated data). We report the results on this par- 
ticular variable because the gas compositions are not mea- 
sured by any sensors and are therefore a potential challenge 
to our algorithm. The error bars represent the uncertainty 
of the estimates as plus and minus two standard deviations 
(for particle filtering we computed the standard deviation 
induced by the weighted samples). 

Although under our setup particle filtering was slower 
than our algorithm by a factor of 750, as Fig. 4 demon- 
strates, the estimates of particle filtering are not as good as 
the estimates of our algorithm. Over the entire sequence the 
average error of our algorithm was 0.000 while the average 
error of particle filtering was 0.01 3. Nevertheless, the more 


dramatic difference is in the estimates of the variance. Of- 
ten, the estimated variance for particle filtering is extremely 
small, even when the estimated value is not very accurate 
(e.g., time steps 72 and 73). In fact, over the entire se- 
quence, according to the estimated distribution of our algo- 
rithm, the correct value of the H? composition was within 
two standard deviations 96% of the time (this is consistent 
with the fact that the probability mass within two standard 
deviations from a Gaussian mean is 95%). In comparison, 
for particle filtering, the true value was within two esti- 
mated standard deviations only 20% of the time. The dif- 
ference was even more apparent when we computed the av- 
erage log-likelihood of the true value, given the two possi- 
ble estimates. For our algorithm the average log-likelihood 
was 3.1 while for particle filtering it was only —5.50 • 1 0 1 1 . 

The reason for this problem is the relatively high dimen- 
sion of the evidence which leads to a very high variance for 
the weights of the samples. Although we generated 10,000 
samples at each time step only a very small number of them 
had a significant effect on the estimate. Over the entire se- 
quence, in 65% of the time steps one sample accounted for 
more than 0.5 of the total probability mass, in 27% one 
sample accounted for more than 0.9 of the mass, and in 
15% one sample accounted for more than 0.99. Obviously 
in cases where one sample completely dominates the rest, 
the estimates of particle filtering are not very reliable and 
in particular the variance estimates can be extremely small 
and misleading. 

Thus, not only is our algorithm faster than particle fil- 
tering with 10,000 samples by a factor of 750, its estimates 
are much more reliable. 

6.2 Results on Real Data 

We next ran a set of experiments on real data. Our data set 
consisted of a long sequence of 13,875 time steps, most of 
it collected while the system was running in steady state. 
We divided our data into a training set, used to estimate 
and tune the model parameters, and a test set on which we 
report our results. 

We conducted a variety of experiments in which we 
compared model predictions with the actual measurements 
recorded by the system under various scenarios: steady 
state and non-steady state, removing sensors, and modify- 
ing the sensor models. In order to make the comparison 
informative, the model predictions for values at time / + 1 
as reported in this section are not adjusted with evidence at 
time /-hi, i.e., they are the predictions based on evidence 
from times 0, 1, . ./. 

Our first experiment, shown in Fig. 5(a), illustrates the 
efficacy of our tracking algorithm during steady-state op- 
eration of the system. In particular, the graph illustrates 
the predicted (thick lines) and measured (thin lines) pres- 
sures, P?, and Pa at point (2) in Fig. 2(a). Observe that 
the predicted value for P?, appears to be consistently lower 
than the measurement. This is the result of the model's 
bias weighting, •) - .07, discussed in Section 4.1, which 
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Figure5: Experimental Results Tracking the RWGS. The .Y axis represents time, and the Y axis the value of the appropriate 
quantity. 


tends to pull the estimates slightly away from the measured 
value. While, in this case, it produces a slightly power re- 
sult, overall, the bias weighting technique does less data 
overfitting and works better in non-steady state sequences. 

Next we experimented with “removing” sensors from 
the system. (This is easily achieved by ignoring selected 
sensor evidence when running the tracking algorithm.) 
Sensor removal can be used to evaluate the robustness of 
the algorithm as well as to determine the importance of a 
sensor for monitoring the system. In Fig. 5(b), we show 
the flow of gas from the compressor at point (11). The two 
overlaid lines are our estimates of this flow value — one 
with all of the sensors, and the other with sensors ft <» and 
R] n (located at ( 13)) removed. In contrast, when flow sen- 
sor R* (located at (16)) is removed, the predicted flow rate 
quickly strays. These results indicate that, at least for this 
sequence, /?s is a more valuable sensor than /??. and R i n. 

We also tested the effects of changing the liquid level 
(LL) sensor noise parameter 4 on our prediction of the gas 
flow ffo at (13). Recall from Section 4.1 that to correctly 
model a sensor we introduced both some Gaussian noise on 
the sensor and a hidden bias variable. We tried both a vari- 

4 The liquid level sensor is very noisy, as splashing and bub- 
bling from the dissolved CO? and from drops splashing from the 
condenser hit the sensor rod and create considerable noise in the 
sensor reading. 


ance value of 0.01 , which we estimated using “reasonable” 
prior knowledge, and a variance value of 4 which was fit to 
the data. Fig. 5(c) shows the effect of the variance of the LL 
sensor for the water tank at (5). With the fitted variance, the 
algorithm tracked quite well. In contrast, with the smaller 
variance, the performance was poor and erratic, following 
the fluctuations in the LL measurements. 

The utility of the bias variables is shown in Fig. 5(d). 
The upper line is a prediction of the flow rate, made using 
a version of the model that contained no bias variables for 
the flow sensors at (10), (13) and (16). The middle line 
corresponds to the model with the bias variables present, 
but shows the prediction for the true (unbiased) flow (i.e., 
the sensor prediction minus the bias). When we explicitly 
modeled the sensor bias, our (unbiased) predictions of the 
true system state better matched the measurements, an in- 
dication of a better estimate of the system state. 

Finally, we tested the ability of the model to track non- 
steady-state behavior — in particular, the behavior of the 
system when the CO? supply is turned off during the shut- 
down process. Unfortunately, we only had one data set con- 
taining this transition, and thus we expect our parameters 
are still not tuned optimally. In addition, having only one 
such transition in our data, we report results on the same 
data that was used for training. 

Fig. 5(e) shows a comparison between the predicted and 







measured output from pressure sensors P?> and Pa, for two 
versions of the model. The first set of predictions, shown in 
solid lines, was calculated using our best estimates of the 
empirical parameters, including the membrane area (calcu- 
lated from other parts of the data set) of 27. 1 . The second 
set of predictions, shown in dashed lines, was calculated us- 
ing an earlier estimate of the membrane area of 3 1 .6. While 
in the steady-state prior to timestep 220, the two predic- 
tions are equivalent as the differences were absorbed into 
the bias errors, in the transient part, the model with inaccu- 
rate parameters underpredicts the initial drop in pressure, 
and retains this error throughout the rest of the sequence. 

Fig. 5(f) presents the predictions of the correct model 
for the flows at ( 16) and fti? ( 10), over a longer period 
of time. Initially, when the CO? supply was cut off, the 
flows dropped; however, gradually the CO and CO? in the 
system were vented and the only remaining gas was H?. As 
the membrane presented less resistance to H? the flow rates 
started to go up. The model tracked this complex behavior 
surprisingly well. 


7 Conclusions and Future Work 

In this paper we address the problem of monitoring a large 
complex physical system — NASA’s Reverse Water Gas 
Shift system — perhaps the largest and most complex hy- 
brid DBN developed to date. This paper makes contri- 
butions both to the modeling and the monitoring of com- 
plex nonlinear systems. On the modeling side, we have 
shown how to model physical systems whose effects man- 
ifest themselves at dramatically different time scales, and 
that involve biased sensors, where the bias is state depen- 
dent and varies over time. On the monitoring side, we have 
presented a general framework for approximating nonlin- 
ear behavior using integration methods that extend the Un- 
scented Filter, improving the accuracy of the approxima- 
tion with minimal additional computation. Experimen- 
tal results indicate that this approach is much faster and 
more reliable than particle filtering. More generally, we 
have demonstrated the feasibility of hybrid DBNs for mon- 
itoring a complex real-world physical system such as the 
RWGS using real data. 

There are several interesting directions for future work. 
The tracking algorithms presented in this paper assume a 
known mode of operation, i.e., all the discrete variables are 
observed. Our long-term goal is to diagnose the RWGS 
when components fail. In order to track both the discrete 
and continuous state, we intend to combine the results pre- 
sented in this paper with algorithms that handle hidden dis- 
crete events such as Rao-Blackwellized Particle Filtering 
(RBPF) [Doucet et al, 2000] or the algorithms presented 
in [Lemer and Parr, 2001; Lemer et al, 2000]. The speed 
of our algorithm (taking just 20ms to generate a Gaussian 
over all the state variables) is a promising indication that 
we can use these techniques for real-time fault diagnosis. 
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Abstract 

Just as actions can have indirect effects on die state of the 
world, so too can sensing actions have indirect effects on 
an agent's state of knowledge. In this paper, we investigate 
“what sensing actions tell us", i.e., what an agent comes to 
know indirectly from the outcome of a sensing action, given 
knowledge of its actions and state constraints that hold in the 
world. To this end, we propose a formalization of the no- 
tion of testing within a dialect of the situation calculus that 
includes knowledge and sensing actions. Realizing this for- 
malization requires addressing the ramification problem for 
sensing actions. We formalize simple tests as sensing ac- 
tions. Complex tests are expressed in the logic programming 
language Golog. We examine what it means to perform a 
test, and how the outcome of a test afreets an agent's state of 
knowledge. Finally, we propose automated reasoning tech- 
niques for test generation and complex-test verification, un- 
der certain restrictions. The work presented in this paper is 
relevant to a number of application domains including diag- 
nostic problem solving, natural language understanding, plan 
recognition, and active vision. 


Introduction 

Agents equipped with perceptual capabilities must operate 
in a world that is only partially observable. To determine 
properties of die world that are not directly observable, an 
agent must use its knowledge of the relationship between 
objects in die world, and its limited perceptual capabili- 
ties to infer such unobservable properties. For example, if 
an agent performs a sense action and observes that there is 
steam coming out of an electric kettle, then the direct effect 
of that sensing action is that the agent knows there is steam 
coming out of die ketde. With appropriate knowledge of the 
functioning of kettles, the agent should also know that the 
electrical outiet has power, that the kettle is functioning, and 
that there is hot liquid inside the kettle - all as indirect ef- 
fects of the sensing action. Similarly, if die agent wishes to 
know whether there is power at an electrical outiet, but can- 
not directly sense this property of the world, die agent may 
potentially acquire this knowledge by attempting to boil wa- 
ter in a ketde plugged into this outlet. 

Copyright © 2000, American Association for Artificial Intelli- 
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Such a sequence of actions constitutes a test. If steam is 
observed, then the agent knows that there is power at the out- 
let; however if steam is not observed, the agent may or may 
not know that there is no power at the electrical outiet. The 
knowledge the agent acquires from the test will depend on 
whether the agent knows that the ketde is functioning. Thus, 
this particular test is only guaranteed to provide knowledge 
about the existence of power at the electrical oudet under 
one test outcome. 

While researchers have extended theories of action to 
include the notion of sensing or knowledge-producing ac- 
tions (e.g., (Scherl Sc Levesque 1993; Baral Sc Tran 1998; 
Golden Sc Weld 1996; Funge 1998)) and have charac- 
terized the effect of sensing actions on an agent's state 
of knowledge, and even how to plan (e.g., (Stone 1998; 
Golden & Weld 1996)) and to project (e.g., (De Giacomo 
St Levesque 1999b)) in certain cases, with sensing actions, 
they have not addressed the problem of how to reason in 
a partially observable environment 1 . More generally, they 
have not examined the problem of how sensing actions 
can be coupled with knowledge of the relationship between 
objects in the world to gain further knowledge, and how 
both sensing actions, and world-altering actions change an 
agents state of knowledge in the presence of such world 
knowledge. Further, they have not examined the problem 
of how to select sensing actions to acquire knowledge of 
some property of the world that is not directly observable. 
Perhaps die closest research is that of (Shanahan 1996b; 
1996a) who investigates the assimilation of sensing results 
for a mobile robot in a framework based on the event cal- 
culus, (Mcllraith 1997) who assimilates observations into 
situation calculus device models to perform dynamical di- 
agnosis, or (Baral, McDraith, & Tran 2000) who do likewise 
in the language C. 

In this paper, we examine these issues in a dialect of the 
situation calculus that has been extended with knowledge- 
producing actions 2 (Scherl Sc Levesque 1993), but which 
does not include state constraints. Following (Mcllraith 
2000), we add state constraints to this language in order to 


1 Partially-Observable Markov Decision Processes (POMDPs) 
address this class of problems within a different formalism, but 
they do not address the testing issues we examine here. 

2 Henceforth referred to simply as sensing actions. 



model the relationship between objects in the world, adopt- 
ing the associated solution to the ramification problem for 
world-altering actions. We show that this solution extends 
to solve the ramification problem in the presence of sens- 
ing actions. Next, we define the notion of a test — how to 
design them and what knowledge can be drawn from their 
outcomes. In the formalization, simple tests comprise a set 
of initial conditions and a primitive sensing action. Complex 
tests are expressed as complex actions in the logic program- 
ming language Golog. We examine what it means to per- 
form a test, and how the outcome of a test affects an agent’s 
state of knowledge. Additionally, we examine the issue of 
selecting tests to confirm, refute, or discriminate a space of 
hypotheses. 

Finally, we investigate the automation of reasoning about 
tests. We show that regression may be used to verify ob- 
jective achievement for complex tests written in a subset of 
Golog. Further restrictions on the form of die complex tests 
allows die same regression operators to serve as the basis 
for a simple regression-style planner that generates tests to 
increase an agent’s knowledge with respect to a space of hy- 
potheses. 

Situation Calculus 

The situation calculus language we use, following (Reiter 
2000), is a first-order language for representing dynamically 
changing worlds in which all of the changes are the direct 
result of named actions performed by some agent, or the in- 
direct result of state constraints. Situations are sequences 
of actions, evolving horn an initial distinguished situation, 
designated by the constant So- If a is an action and s a sit- 
uation, the result of performing a in s is the situation rep- 
resented by the function do(a , s). Functions and relations 
whose truth values vary from situation to situation, called 
fluents , are denoted by a predicate symbol taking a situation 
term as the last argument. Note that for the purposes of this 
paper, we assume that our theory contains no functional flu- 
ents. Finally, Posh {a. s) is a distinguished fluent expressing 
that action a is possible to perform in situation s. A situation 
calculus theory V comprises the following sets of axioms: 

• foundational axioms of the situation calculus, E, 

• successor state axioms, V$s, 

• action precondition axioms, V ap , 

• axioms describing the initial situation, T>s n , 

• unique names for actions, V un a* 

• domain closure axioms for actions, 

Successor state axioms, originally proposed by (Reiter 
1991) to address the flame problem and extended by (e.g., 
(Lin & Reiter 1994; Mcllraith 2000)) to address the ramifi- 
cation problem, are created by making a causal interpreta- 
tion of die ramification constraints and a causal complete- 
ness assumption and compiling effect axioms of the fonn 3 : 

Pons {a. s) Ayr (*? « s *) 3 F(x, do(a , ,s) ) (1 ) 

Poss(a, s) A 7 p (*-, «» *) 3 “'Ffl. do(a. .s)) . (2) 

3 Notational convention: all formulae are universally quantified 
with maximum scope unless otherwise noted. 


and ramification (state) constraints of the form: 

v+{x,s)OF{x,s) (3) 

D ->F(x,s\ (4) 

into Intermediate Successor State Axioms of the form: 

Poss(a.s) 3 \F t (x. do{a. s)) = <bp.] where, (5) 

<&r, = 7f; <*, «> •’) V i'f, (••?, ■•>)) 

V(F(?,s) 

A -'(7F, (*> «, •’) v i.y. (.if, i lot a, . 1 )))), (6) 

I.e., if an action is possible is situation s, then it implies that 
die fluent is true in do(a , .*) iff an action made it true -or- 
a state constraint made it true -or- it was already true and 
neither an action nor a state constraint made it false. 

Such intermediate successor state axioms provide a com- 
pact representation of a solution to the ramification problem 
for a common class of state constraints. (Mcllraith 2000) 
shows that for what are essentially acyclic causal ramifi- 
cation constraints, repeated regression rewriting (e.g., (Re- 
iter 1991)) of$p.,7v* [$/?.] = repeatedly rewrites the 
r amific ation constraints that are relativized to do(a , .s) in (6) 
above, and is guaranteed to terminate in a formula whose 
fluents are relativized to situation s rather than do(a.s). 
Both die intermediate and the less compact (final) succes- 
sor state axioms which result from die regression provide 
closed-form solutions to the flame and ramification problem 
for the designated class of state constraints. 

To illustrate sensing and testing in partially observ- 
able environments, we present a partial axiomatization of 
a car repair domain, derived from The Complete Idiot's 
Guide to Trouble-Free Car Care (Ramsey 1999). Our do- 
main includes world-altering actions such as tumj>n(x) and 
turn -of fix), where x is radio or lights. These have the 
effect that the radio or lights are on/off in die resulting situ- 
ation. Actions turn(key) and released key) have the effect 
that the ignition is begin turned (turning Jgn) 9 or not, in 
the resulting situation. These actions are defined in terms 
of effect axioms and are combined with the following self- 
explanatory state constraints to produce successor state ax- 
ioms. For notational convenience we abbreviate: transmis- 
sion - travs, interlock - intrlk , solenoid - solnd, engine - 
engn, battery - haft; ignition system - ignsys, start system - 


strtsys. 

cmptyigafiJank.fi) 7 ->) itartableis) (7) 

abiintrlk.fi) 3 ->fitartableis) (8) 

ubibatt.fi) 3 -miartableifi) (9) 

ub(fiofntL fi) 3 -ifiiarlnhteifi) (10) 

abifitarier.fi) 3 ^fifar tablets) (11) 

au/oUrnn*) A ~>ingear{iranfi. *) 3 ab{inlrik. fi) (12) 
manual(l rans) A -aleprefified{clutch, fi) 3 abiintrlk.fi) (13) 
turning Jgnifi) A ab{batt , *) 3 ->noific(engn. fi) (14) 
turning Jgnis) A einptyigafiJunk.fi) 3 ->noifie{engn. *) (15) 



f timing Jgn{s) A ~^nti[aohid. #) D noijn y (*olnd. *) (16) 

abibati. j») A <m( radio . *) 3 -inoi#e( radio , *) (17) 

n/»( radio, 3 (18) 

->ab(balLn) Aon{fighU. a) 3 emit &(tight < s ) 0$) 


Space precludes listing all the successor state axioms. There 
is one (intermediate) successor state axiom for each fluent. 
E.g., axioms (7H11) compile into intermediate successor 
state axiom (20): 

Pos$(a, s) D [startaUe(do(a t s)) s 

-^cmptyiga sJta nk. do(a , s)) A -tabiintrlk. do(a. s)) 

A -i abibatt , </o(tt, ,«*)) A ^<zfc(s0f7*/.do(a, -s)) 

A -*ab(startcr. do(a. ,s))] (20) 

As described in (Mcllraith 2000), the axioms describing 
the initial situation, So contain what is known of the initial 
situation as well as the ramification constraints of the form 
of (3) and (4), relativized to Sq. 

Knowledge and the Ramification Problem 

In (Scherl & Levesque 1993), the situation calculus lan- 
guage without state constraints was extended to incor- 
porate both knowledge and sensing actions. World- 
altering actions change the state of the world, sensing ac- 
tions have no effect on the state of the world but rather 
change the agent’s state of knowledge. In our exam- 
ple, sensing actions include check -f uel* check ^ar Mart, 
check. radio jnoisc etc., which have the effect of the agent 
knowing empty (gas Jank, do(a % *)), st actable (<lo(a, s)) t and 
noise ( radio , do(a, s)). 

The notation Knows(o,s) (read as o is known in 
situation s) y where <i> arbitrary formula, is an abbre- 
viation for a formula that uses K. For example 
Knows(o7?.( block ] , block?)* s) abbreviates: 

V.s' K(fi , s) 3) onlbloek) , block ? , d). 

The notation Kwhether(<y, ,s) is an abbreviation for a for- 
mula indicating that die truth value of o is known. 

Kvrlietlier(^, s) = r Knowsia, s) V Knowsf-^. 

Following the notation of (Levesque 1996), each sense ac- 
tion a has a sensed fluent , SF(a. s) associated with it, and 
for each such a, T> entails a sensed fluent axiom: 

SF(a,s) = vi(s), (21) 

which says that performing the sense action a tells the agent 
whether die formula p(s) is true or false. Thus, T> 
Kwhether ( {$?, do(s, s)) where a is an action with a sensed 
fluent equivalent to o. 

For the sense action check-f uel the sensed fluent axiom 
is: 

SF (check .fuel, s) = empty (gas Ja nk. s) (22) 

which tells us whether or not die gas tank is empty. For 
world-altering actions, V entails SF(a , .s) = True. 

In (Scherl & Levesque 1993), a successor state axiom for 
the K fluent is developed. Its form is as follows: 

Successor State Axiom for K 
Poss(a.s) D [K(s" .do(a, s)) = 

3 s'. Pass (a. s') A K(s\s) A (s' — do(a. s')) A 
[$F(o,y)sSF(a,*)] (23) 


which says that after doing action a in situation .s, the agent 
thinks it could be in a situation s n iff .s" = do(a. s f ) and s' 
is a situation that was accessible from .% and where s and *' 
agreed on the truth value of SF(n, s ), e.g., the truth value 
of empty (gas Jank). Thus, for all situations <lo(a. s), the K 
relation will be completely determined by the K relation at $ 
and the action a. This extends Reiter’s solution to the frame 
problem (without ramifications and without knowledge) to 
the case of the situation calculus with sensing actions. 

Proposition 1 In the situation calculus theory described 
above, the agent knows the successor state axioms and the 
ramification constraints. 

This follows from the fact that die successor state axioms 
are universally quantified over all situations, and the rami- 
fication constraints explicitly hold in So and are entailed in 
all successor situations, by the successor state axioms. 

Theorem 1 (Correctness of Solution) The proposed solu- 
tion to the frame and ramification problems for world- 
altering and sensing actions ensures that knowledge only 
changes as appropriate , as defined by Theorems 1, 2, 
3 (Scherl & Levesque 1993). Furthermore , the agent knows 
the indirect effects of its sensing actions. 

Thus, the successor state axioms for world-altering and sens- 
ing actions, together address the frame and ramification 
problems. 

Testing 

The purpose of a test is to attempt to determine the truth 
value of certain properties of the world, that may or may 
not be directly observable. A test is often performed with 
respect to a set of hypotheses, with the objective of elimi- 
nating as many hypotheses as possible from the set of hy- 
potheses being entertained Testing has been studied ex- 
tensively for the specific problem of IC circuit testing, but 
there is little work on testing for rich dynamical systems 
such as the ones we examine here. The notion of a static 
test was briefly discussed in (Moore 1985, litmus example), 
and further developed for static systems in (Mcllraith 1994; 
Mcllraith & Reiter 1992). We build directly upon the work 
in (Mcllraith 1994) with the objective of developing a for- 
mal theory of testing for dynamical systems. 

Informally, a simple test comprises a set of initial con- 
ditions that may be established by the agent, together with 
the specification of a primitive sensing action, which deter- 
mines what the agent will directly come to know as the result 
of the test. In our car repair domain, we can test the battery 
by checking the radio for noise. The initial conditions for 
such a test might be on ( radio, s). Then we can perform the 
sensing action cheek j'adiojnois e to see whether the radio is 
emitting noise. Note that the precondition for petforming 
the action check jrtMliojnoi sc , Poss(che ckjradiomoi sc, s) ~ 
i n side (car, s), is different from the initial conditions of the 
test. Both must hold and must be consistent with the theory 
and with the current hypotheses being entertained, in order 
to execute the test. 

We distinguish between two types of tests, truth tests 
which tell us whether the properties being sensed are true in 


the physical world, and junctional tests , which tell us what 
values of the properties are true in the physical world. For 
the purposes of this paper, we restrict our attention to truth 
tests, and our sensing actions to so-called binary sense ac- 
tions which establish the truth or falsity of a sensed formula. 

Definition 1 (Simple Test) 

A simple test is a pair, (I, a), where I, the initial conditions , 
is a conjunction of literals, and a is a binary sense action 
whose sensed formula contains no free variables. 

(oniradio, s), check j'adiojnoisc) is an example of a simple 
test, following the discussion above. We now define the no- 
tion of a test for a particular hypothesis space, represented 
by the set HYP. We restrict the hypotheses, His) £ HYP 
to be conjunctions of fluents whose non-situation terms are 
constants, and whose situation term is a situation variable s. 
In our car repair domain, an example hypothesis space might 
be {ab(batt, s ) , ah(soliul, ,s). cmptyigasJank, .s)}. 

Definition 2 (Test for Hypothesis Space HYP) 

A test (/, a) is a test for hypothesis space HYP in situation 
s iffV AlAPofss{a , s)aH(s) is satisfiable for every H(s) £ 
HYP . 

That is, the state the world must be in to execute 
the sensing action must be satisfiable, under the as- 
sumption that any one of the hypotheses in die hypoth- 
esis space could be true. Consider that V entails the 
safety constraint -exploswn(s) and the axiom sparks(s) A 
gasJcak(s) 3 explosion(s\ and that our hypothesis space 
is {gasJeak(s),ab(sparkj)lug, s)}. A reasonable test for 
ab( spark-plug, s) is to try to create sparks at die plug. Unfor- 
tunately such a test would cause an explosion in die presence 
of a gas leak. The satisfiability check above precludes such 
a test. 

Definition 3 (Confirmation, Refutation) 

The outcome a of the test {La) confirms H(s) 6 HYP 
iff V A 1 A Poss(a , s) A H(s) is satisfiable and T> A I A 
Poss(a. s) \= Knows(i7 D a.s). a refutes H(s) iffV A 
I A Poss (a. s)aH(s) is satisfiable and V A 7 A Post (a. s) |= 
Knows (H D -««, s). 

If the outcome of test (on (radio, s), check -radio jnoisc) is 
noise {radio, do(a. *)), then our test refutes die hypothesis 
abibatt . ,s), following Axiom (17), and we can eliminate 
abilMitt, s) from our hypothesis space, HYP . 

Observe that a test outcome that refutes an hypothesis 
H(s) allows us to eliminate it from HYP. Unfortunately, a 
test outcome that confirms an hypothesis is generally of no 
deterministic value, resulting in no reduction in the space of 
hypotheses. As we will see in a section to follow, there are 
exceptions that depend on the criteria by which die hypoth- 
esis space is defined. 

In the sections to follow we use these basic definitions 
to define discriminating tests and relevant tests. These tests 
are distinguished by the effect their outcome will have on a 
general space of hypotheses. 

Discriminating Tests 

Notice that in our example above, if we had observed 
^iwisciradio, do(a, s)), then by die definition, this would 


have confirmed the hypothesis abibatt, s), but it would have 
been of little value in discriminating our hypothesis space. 
All hypotheses remain in contention. Discriminating tests 
are those tests (7, a) that are guaranteed to discriminate an 
hypothesis space 77YP, i.e., which will refute at least one 
hypothesis in HYP, regardless of the test outcome. 
Definition 4 (Discriminating Tests) 

A test (L a) is a discriminating test for the hypothesis space 
HYP iffV A 7 A Poss(a , .s) A H(s) is satisfiable for all 
H(s) e HYP, and there exists Hi(s), Hj(s) € HYP 
such that the outcome a of test (7, a) refutes either Hj{s) 
or Hj (s), no matter what that outcome might be. 

Proposition 2 

After we perform a discriminating test, (L a), 
Knows(-i77j , $), for some Hj (*) € HYP. 

In general, we would like a discriminating test to reftrte 
half of the hypotheses in the hypothesis space, regardless of 
the test outcome. By definition, a discriminating test must 
refute at least one hypothesis in the hypothesis space. 

Definition 5 (Minimal Discriminating Tests) 

A discriminating test (7. a) for the hypothesis space HYP 
is minimal iff for no proper subconjunct V of I is (7 / .a) a 
discriminating test for HYP. 

Minimal discriminating tests preclude unnecessary initial 
conditions for a test. 

In some cases, we are interested in identifying a test that 
will establish the truth or falsity of a particular hypothesis. 
An individual discriminating test does precisely this. 

Definition 6 (Individual Discriminating Tests) 

A test (7, a) is an individual discriminating test for the hy- 
potheses Hi(s) and ->Hi(s) € HYP iffVAlAPoss(a, s)A 
H{s) is satisfiable for all H(s) £ HYP and the outcome a 
of test (I. a) refutes either Hi(s) or -<Hi(s), no matter what 
that outcome might be. 

Proposition 3 

After we perform an individual discriminating test (La), 
Kwhether(ifi. s) for some Hi € HYP. 

The test ({}, check -f uel) is such a test. The out- 
come will be one of -*empty(gasJank, do(a, .s)) or 
cmptyigasJank, do(a,s)). Thus, as the result of per- 
forming check-fuel in the physical world, the agent 
KwhetixtTicmpt y(ga*Jank, s)). 

We can similarly define the notion of a minimal individual 
discriminating test, and a minimal relevant test, below. 

Relevant Tests 

In the majority of cases we will not be so fortunate as to 
have discriminating tests. Relevant tests are those tests 
(7. a) that have the potential to discriminate an hypoth- 
esis space HYP , but which cannot be guaranteed to do 
so. Given a particular outcome cv, a relevant test may re- 
fute a subset of the hypotheses in the hypothesis space 
HYP, but may not refute any hypotheses if is ob- 
served Since we can’t guarantee the outcome of a test, 
these tests are not guaranteed to discriminate an hypothe- 
sis space. (oniradio, .<*). check j'OAliojnoise ) is an example of 
such a test. 



Definition 7 (Relevant Tests) 

A test {La) is a relevant test for the hypothesis space 
HYP iffV A / A Poss(a , s) A iff*) is satisfiable for all 
H {.<t)inHY P, and the outcome a of test (L a) either con- 
firms a subset of the hypotheses in HYP or refutes a subset. 

By definition, a relevant test confirms or refutes at least 
one hypothesis in HYP, and it follows that every discrimi- 
nating test is a relevant test. 

In addition to discriminating and relevant tests, there is 
a third class of tests. Constraining tests do not refute an 
hypothesis, regardless of the outcome, but they do provide 
further knowledge drat is relevant to the hypothesis space 
and which the agent can exploit in combination with other 
tests. We discuss this notion in a longer paper. 

Testing Hypotheses 

In the previous section we observed that a test outcome that 
refutes an hypothesis H{s) e HYP allows us to eliminate 
it from HYP, but that in general an outcome that confirms 
H{s) has no value in reducing die hypothesis space. In this 
section, following (Mcllraith 1994), we show that when the 
hypothesis space is determined using a consistency-based 
criterion this is indeed true, but when the hypothesis space is 
defined abductively, confirming test outcomes serve to elim- 
inate those hypotheses that are not confirmed, i.e., that do 
not explain, the test outcome. 

Definition 8 (Consistency-Based Hypothesis Space) 

A consistency-based hypothesis for V and outcome o of 
the test {I, a) is any H (a) € HYP such that V A / A 
Poss(a , s) A H(s) Act is satisfiable . 

Proposition 4 (Eliminating C-B Hypotheses) 

The outcome a of a test {L a) eliminates those consistency- 
based hypotheses , H{s) € HYP that are refuted by test 
outcome rv. 

Definition 9 (Abductive Hypothesis Space) 

An abductive hypothesis for V and outcome a of the test 
(/. a) is any H{s) € HYP such that V A / A Poss{a , .$) A 
H(s) is satisfiable, and V A I A Poss{a , .s) A H(s) f= a. 

Proposition 5 (Eliminating Abductive Hypotheses) 

The outcome a of a test (L a) eliminates those abductive 
hypotheses, H(s) e HYP that are not confirmed by test 
outcome a. 

Thus, in the case of abductive hypotheses, unlike 
consistency-based hypotheses, both confirming and refuting 
test outcomes have die potential to eliminate hypotheses. 

Proposition 6 (Efficacy of Tests) 

Any outcome o of a relevant test (I.a) can eliminate abduc- 
tive hypotheses, whereas only a refuting outcome can elimi- 
nate consistency-based hypotheses . Discriminatory test out- 
comes, by definition, can eliminate either consistency-based 
or abductive hypotheses, regardless of the outcome . 

Complex Tests 

In the previous section, we defined the notion of a simple 
test (/. a), and characterized the circumstances under which 


the outcome of such a test would discriminate an hypoth- 
esis space. Indeed, to discriminate an hypothesis space, we 
may need a sequence of simple tests, interleaved with world- 
altering actions in order to achieve the initial conditions for 
a test. Likewise, the selection and sequencing of sensing 
and world-altering actions may be conditioned on the out- 
come of previous sensing actions. In the section to follow, 
we examine the problem of generating tests using regres- 
sion. As we will see, generating tests, especially tests that 
involve sequences of sensing and world-altering actions is 
hard. In many instances, we need not resort to computation. 
The domain axiomatizer can articulate procedures for testing 
aspects of a system, just as the author of The Idiot's Guide 
has done in the domain of car repair. The logic programming 
language, Golog (alGOl in LOGic) (Levesque et al 1997) 
provides a compelling language for specifying such tests, as 
we describe briefly here. 

Only a sketch of Golog is given here. See (Levesque et al. 
1997) for a full discussion of the language and also a Prolog 
interpreter. Golog provides a set of extralogical constructs 
(such as action sequencing, if-then-else, while loops) for as- 
sembling primitive actions, defined in the situation calculus, 
into macros that can be viewed as complex actions. The 
macros are defined through the predicate Do{6, s, .s') where 
t) is a complex action expression. Do{ <\ s, s') is intended to 
mean that the agent’s doing action <) in situation s leads to 
a (not necessarily unique) situation s'. The inductive defini- 
tion of Do includes the following cases: 

Do(a. s. s') — simple actions 

Do(6?, s. s) — tests (referred to as G-tests in this paper) 
Do( [S i ; S- 2 ] , s, s ) — sequences 
Do( [S] s. s') — nondeterministic choice of actions 

Do{(Th:)S. s. s') — nondetenninistic choice of parameters 
DoiiS <p then S\ else S>, s, s')- conditionals, where we 

restrict o toaG-test 
DofwhUe <p do<L,s.s') — while loops 

Space does not permit giving file lull expansion far each 
of the constructs, but they can be found in (Levesque et al. 
1997). The only change here is that the definition of the G- 
test construct (including the implicit G-test in file condition 
construct) must expand into a G-test involving knowledge 4 . 

The following is a partial example of a complex test writ- 
ten in Golog, and derived from (Ramsey 1999). This par- 
ticular procedure is designed to help discriminate the space 
of hypotheses generated when a car won’t start, namely 
{ab(intrlk, s), empt y ( gas J ark. ,s). ab(batt, ,s), ab(solrul , s). 
ab(igji .wires,*), ab{ starter,*)}. In a diagnostic application 
such as this one, Golog procedures may also be written to 
combine testing with repair. 

proc CarWontStart 
I f (-i startable) then checkInterlock; 

4 We are taking the simplest approach towards incorporating 
sensing actions into Golog. All actions are on-line. In other words, 
they are executed immediately without any possibility of back- 
tracking. Other options for completely off-line execution (Lake- 
meyer 1999) and a mixture of off-line and on-line execution (De 
Giacomo & Levesque 1999a) have been discussed in the literature. 


if (— . AB(lNTRLK)) then CHECK JGAS_TANK; 
if (-< empty(gas.tank)) then checkBattery; 
if (-> ab(batt)) then checkSolenoid; 
if (-» ab(solnd)) then checkIgn Wires; 
if (-h ab(ign_wires)) then checkStarter; 
if (-1 ab(starter)) then checkEngine 
end if end if end if end if end if end if end if 
endProc 

proc CheckBattery 
turn^on(radio); check_radio_noise; 

if (-. NOISE(RADIO)) 

then turnjon(lights); check_lights 

end if 
endProc 

Observe that complex tests often invoive world-altering 
actions which serve to establish the preconditions and initial 
conditions for embedded simple tests. Also observe that in 
achieving the preconditions or initial conditions for simple 
tests, these actions change the state of die world, including 
potentially changing the space of hypotheses. For exam- 
ple, if a flashlight isn’t emitting light, and one hypothesis 
is that the batteries are dead, a good way to test them is to 
replace them with fresh batteries, and see whether the flash- 
light then works. However, replacing the flashlight batteries 
potentially changes the state of one of the hypotheses. 

In diagnosis domains, such as the ones above, it is of- 
ten desirable to combine fault detection (hypothesis testing) 
with repair and to take actions to eradicate faults as easily as 
to diagnose them (Mcllraith 1997; Baral, Mcllraith, & Tran 
2000). However, in cases where it is desirable not to alter 
the truth status of the hypothesis space, care must be taken 
to design and verify and/or generate tests that maintain des- 
ignated knowledge constraints and world constraints. E.g., 
we don’t want to determine whether the gas tank is empty 
by draining it! 


involve ensuring that for at least one H, Knows(->XL s) 
holds in the final situation, i.e., V/^m r Knows(->/f. s). 

In (Scherl & Levesque 1 993), a form of regression (based 
on the discussion in (Reiter 1991)) is developed for the sit- 
uation calculus with sensing actions. Through the appli- 
cation of regression, reasoning about situations reduces to 
reasoning in the initial situation, So* Given a ground sit- 
uation term (i.e. a term built on So with the function do 
and ground action terms) the problem is to determine 
whether the axiomatization of the domain V entails G(s gr ) 
where G (the intended objective of the procedure) is an arbi- 
trary sentence including knowledge operators. This question 
is reduced to the question of whether or not the axiomatiza- 
tion of the initial situation entails the regression of G(s gr ), 

i. e.,ft(G'( Agr))- Since the result of regression is a formula in 
an ordinary modal logic of knowledge (i.e. a formula with- 
out action terms and where die only situation term is So) an 
ordinary modal theorem proving method may be used to de- 
termine whether or not the regressed formula is entailed by 
the axiomatization of the initial situation, V.%. In our case 
G will be a formula made up of subformulae of the form 
Kwhether(/L a) or Knows (-iff. a), where H is an hy- 
pothesis. 

The regression operator ft is defined relative to a set of 
successor state axioms XL*- The first four parts of the defi- 
nition of the regression operator 5 , ft concern world-altering 
actions and are taken from (Reiter 2000). 

L When W is a non-fluent atom, including equality atoms, and 
atoms with the predicate symbol Poss, or when IT’ is a fluent 
atom or Knows operator, whose situation argument is the situa- 
tion constant So, ft [II I = U\ 

ii. When F is a relational fluent (other than K) atom whose suc- 
cessor state axiom in />«?.$ is 

, Possia , s) D [F(x \ ..... , do(a. a)) = # r] 


Automated Reasoning About Tests 

In the previous section we introduced the notion of a com- 
plex test, demonstrating that such tests could sometimes be 
specified in Golog. In this final technical section we briefly 
examine the use of automated reasoning techniques, and in 
particular the use of regression rewriting, for the purpose 
of verifying certain properties of Golog-specified complex 
tests, and for generating complex tests as conditional plans. 
Our presentation draws upon (Lesperance 1994) and (Re- 
iter 2000). Other related approaches to conditional plan- 
ning include (Rosenschein 1981 ; Manna & Waldinger 1987; 
Lobo 1998). 

Consider the Golog complex test given above to help dis- 
criminate the space of hypotheses generated when a car 
won’t start. To verify that it is an individual discriminat- 
ing test, it is necessary to ensure that for at least one of the 
hypotheses H , Kwhether(fiL a) holds, where a is the sit- 
uation resulting from the execution of the Golog procedure, 
i.e. Do(CarWontStart,5o,a). Thus, we would like to 
be able to entail Vw<=//y r Kwhether(/L a), and in par- 
ticular K whet her(emply(gaAjank), a), for example. A 
verification that the procedure is a discriminating test would 


ft[F(f,, ... 


= if;:.: 


.a.* 


lit Whenever IT ’ is a formula. 


iIT'1 = -ft[IT'] ; 
(Vv)\V] = <Vi.)ft[tr], 
‘(3t OIL,] = (3t,.)ft[iri]. 


iv. Whenever W\ and \V 2 are formulas, 
ft[ILi A II 2] = ftflT'i] A ft|UU 
ftpT’i v ii 2] = ftpr,] v ft.[ir, l 
ftir. Dir 2 ] = ft[ir,]Dft.[ir,] 


Following (Scherl & Levesque 1993), additional steps are 
needed to extend the regression operatorto sensing actions 6 . 
Two definitions are needed for the specification to follow. 
When ^ is an arbitrary sentence and $ a situation term, then 
V?[s] is the sentence that results from adding an extra argu- 
ment to every fluent of ■•+> and inserting a into that argument 

5 Some details are omitted here (e.g, regression of functional 
fluents, and the equality predicate). Also note that the formula to 
be regressed must be regressable. This concept is fully defined in 
(Reiter 2000). 

6 Regression of sensing actions that make known the denotation 
of a term (e.g. an action of reading a number on a piece of paper) 
is not discussed here. 



position. The reverse operation 1 is the result of remov- 
ing the last argument position from all the fluents in p. 

Step v covers the case of regressing a world-altering ac- 
tion through the Knows operator. Step vi covers the cases of 
regressing a sensing action through the Knows operator. In 
the definitions below, a' is a new situation variable. 

v. Whenever a is not a sensing action, 

7? [Know* (IF, <to(o.. a))] = 

Knows ((ft[IF[<fo(a, a')]])~ 1 , a). 

vt Whenever a is a sensing action, where is a formula 
such that Z> entails that V’M 15 equivalent to SF(a,n ), 
77[Knows(TF.</o(a, a))] = 

(Its (a) 3 Knowi(r» t 3 ft[lF[do(a, a')]] ‘.a)) A 
3 Kitows(-»r\ 3 7l[W[d<Aa. V)]]~\.*» 

An additional operator C needs to be defined to handle 
the expansion of me complex actions found in Golog, so 
that we can apply regression 7 . We are only considering a 
subset of Golog programs - those composed of simple ac- 
tions, sequencing, and conditionals. We also add the empty 
action noOp or [j (names for the same operation). Also note 
that 7r a (i\ s) stands for the preconditions of a(:r) as speci- 
fied in the action precondition axiom, V ap , Po*a («(£), s) = 

vtU. CfnoOp, IF, s) = «••(,) 

lx. (J([u(iO; IF, a) = n a (x : a) A C(6, IF, do(a(x), a)) where 
a(.r) is a ground non-sensing simple action term. 

X. F([if <p (.tT) then S j else S - i ], IF, a) = 

K whet her ( o(.?). a) A 
[Knows(^(/) r ?) 3 , IF, a)] A 

[Knows! -v>(£). .$) 3 C(Si, IF, a)] 

We are assuming that the agent is able 8 to execute the Golog 
test procedure. In particular, the programmer (of the test 
procedure) must have ensured that at the point where an 
[if (-) (?) then 6 \ else A- 2 ] statement is encountered, the ex- 
ecuting agent must K whet her (o, a). If not, the procedure 
will Ml. 

In the following theorem (a generalization of Theorem 2 
from (Lesperance 1994), recall Tv* ( p) indicates the repeated 
regression of p until further applications leave the formula 
unchanged. 

Theorem 2 For any Golog procedure 6, consisting of sim- 
ple actions , sequences , and conditionals , and G an arbitrary 
closed regressable formula that may include knowledge op- 
erators: 

T>£= 3s(Do(6.Soi s) A (7(a)) iff 

x>A 0 u x >„„ 0 >= nr ids, g, s*)) 

Theorem 2 shows it may be verified that any Golog testing 
routine (utilizing concatenation and conditionals) achieves 
its intended objective G through the use of regression fol- 
lowed by theorem proving in the initial database. The suc- 
cessor state axioms (!>,,) are only used in the regression 
procedure. This theorem can be extended to likewise verify 
other properties of our Golog procedures. 

7 The C operator introduced here is based on (but generalizes) 
the E operator of (Lesperance 1994). 

8 See (Lesperance et al 2000) for a discussion of ability and 
Golog programs. Related issues are discussed in (Lesperance 1994; 
Lakemeyer 1999). 


We can use the above regression operator as the basis 
for a simple conditional planning algorithm for constructing 
complex tests. Following (Lesperance 1994), we consider 
only normal form conditional plans. These are conditional 
plans in which the condition in a conditional (e.g. the e> in 
[if O (/) then rfi else rfo]) must be a sensed formula. Thus 
we can require that prior to any conditional with the G-test 
o, there must be an action a such that a is a sensing action 
and V (= SF(a, a) = c>(a). This guarantees that the pro- 
gram executing the test will always Kwhether(o, a) when 
a conditional is encountered. For any complex test (that is 
executable) consisting only of concatenation and condition- 
als, there must be an equivalent test in this normal form. 

For i = 1, 2, 3 . . . , we can define the sentences Ti as: 

r 0 =G(») 
r /'= r 

3a([3.r (a = .1, (j*) A , (.r) V ... 

V 3? (a ;= A„{£) A -4 fJ (iO)] 

a ft(r,_i(,/o(«,A))}) v 

3a([3.r (a = - A x A * (x) A (SF{a. a) ~ A 

71(0 \ (J?, do(a. a)) 3 I\_ i (do(a. a))) A 
72(^o,(.r./fo(a.*)) 3 h-i (do(a. a)))] 

A... A 

3a([3£(a - . At r A * t (/) A (SF(a, a) = <$,*(*)) A 

7 (*, do(a, a)) 3 r,_ 1 (do(a, s))) A 
7l(^o m (.f. <io(a , a)) 3 IV 1 (f tola, a)))] 

Each T, is true if there is a plan of length i starting in a 
and leading to a state satisfying G (Reiter 1995; Lesperance 
1994). The following theorem (essentially Theorem 3 of 
(Lesperance 1994)) establishes the soundness and complete- 
ness of the regression-based test planning method. 

Theorem 3 For Golog procedure 6 in normal form and G, 
an arbitrary closed regressable formula that may include 
knowledge operators: 

V [= 3s(Do(H. So, a) A (7(a)) iff for some n 
T>s n UT>„„ a (= To(5 0 ) V ... V r„(S 0 ) 

This regression-based finite horizon method of generating 
and evaluating all normal form conditional plans of greater 
and greater size is certainly not designed for efficiency, but 
the results can serve as the foundation for building more 
efficient regression-based complex-test planning methods, 
much as similar results have served as the foundation for 
relatively more efficient regression based planning methods 
(McDermott 1991 ; Lesperance 1994; Rosenschein 1981). In 
future work we will evaluate the extension of current state 
of the ait planning techniques based on SAT and Graphplan, 
to address the planning problems raised in this paper (Weld 
1999). 

Summary 

In this paper we presented results towards a formal theory 
of testing for dynamical systems, specified in the language 
of the situation calculus. Our first contribution was to ad- 
dress the ramification problem for sensing actions. We then 
defined the notion of a test, examining how a test can be 
designed and how die outcome of different types of tests af- 
fect an agent’s state of knowledge. The realization of many 



tests in the world requires a complex sequencing of world- 
altering and sensing actions, whose selection and ordering is 
conditioned upon the outcome of previous sensing actions. 
We proposed specifying such complex tests in the logic pro- 
gramming language Golog. We then demonstrated that re- 
gression could be used both to verify the desired objective 
of such complex tests, and to generate tests as conditional 
plans under certain restrictions. 

Sensing is integral to die operation of most autonomous 
agents. The notion of complex and simple tests introduced 
here extends the body of theoretical work on sensing in dy- 
namical systems, and has practical relevance for building 
agents for diagnostic problem solving, plan understanding, 
or simply for mobile cognitive agents that need to interact in 
complex environments with limited sensing. 
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Abstract 

In this paper we examine the problem of monitoring and di- 
agnosing noisy complex dynamical systems that are modeled 
as hybrid systems - models of continuous behavior, inter- 
leaved by discrete transitions. In particular, we examine con- 
tinuous systems with embedded supervisory controllers that 
experience abrupt, partial or full failure of component de- 
vices. Building on our previous work in this area (MBCG99; 
MBCGOO), our specific focus in this paper is on the mathe- 
matical formulation of the hybrid monitoring and diagnosis 
task as a Bayesian model tracking and selection problem, and 
provision of a suitable tracking algorithm. The nonlinear dy- 
namics of many hybrid systems present challenges to prob- 
abilistic tracking. Further, probabilistic tracking of a system 
for the purposes of diagnosis is problematic because the mod- 
els of die system corresponding to failure modes are numer- 
ous and generally very unlikely. To focus tracking on these 
unlikely models and to reduce the number of potential mod- 
els under consideration, we exploit logic-based techniques for 
qualitative model-based diagnosis to conjecture a limited ini- 
tial set of consistent candidate models. In this paper we dis- 
cuss alternative tracking techniques that are relevant to dif- 
ferent classes of hybrid systems, focusing specifically on a 
method for tracking multiple models of nonlinear behavior 
simultaneously using factored sampling and conditional den- 
sity propagation To illustrate and motivate the approach de- 
scribed in this paper we examine the problem of monitoring 
and diagnosing NASA’s Sprint AERCam, a small spherical 
robotic camera unit with 12 thrusters that enable both linear 
and rotational motion. 

Introduction 

We have been conducting an ongoing project to investigate 
how to diagnose hybrid systems - complex dynamical sys- 
tems whose behavior is modeled as a hybrid system. Follow- 
ing die description in (MBCG99; MBCGOO), hybrid mod- 
els comprise both discrete and continuous behavior. They 
are typically represented as a sequence of piecewise con- 
tinuous behaviors interleaved with discrete transitions (e.g., 
(Bra95)). Each period of continuous behavior represents a 
so-called mode of the system. For example, in die case of 
NASA’s Sprint AERCam, a spherical airborne robot cam- 
era unit, modes might include translate X-axis, rotate X- 
axis, translate X-axis, etc. (AG98). In the case of an Airbus 
fly-by-wire system, modes might include take-off, landing , 
climbing , and cruise. Mode transitions generally result in 


changes to the set of equations governing the continuous be- 
havior of the system, as well as to die state vector that initial- 
izes tiiat behavior in the new mode. Discrete transitions that 
dictate such mode switching are modeled by finite state au- 
tomata, temporal logics, switching functions, or some other 
transition system, while continuous behavior within a mode 
is modeled by, e.g., ordinary differential equations (ODEs), 
difference equations, or differential and algebraic equations 
(DAEs). For the purposes of this paper, we restrict our at- 
tention to discrete-time estimation for the class of systems 
whose hybrid models contain no autonomous jumps. I.e., 
all nominal transitions between system modes are induced 
by a controller action; none are induced by the system state 
and mode (Bra95). 

In (MBCG99) we presented the hybrid diagnosis prob- 
lem: 

Given a hybrid model of system behavior, a history of 
executed controller actions, a history of observations, 
including observations of aberrant behavior relative to 
the model, isolate the fault that is the cause for the 
aberrant behavior. 

Our task was to perform diagnosis online in conjunction 
with the continued operation of the system. Hence, we 
divided our diagnosis task into two stages, initial conjec- 
turing of candidate diagnoses and subsequent refinement 
and tracking to select the most likely diagnoses. We cast 
die diagnosis problem as die problem of finding a model 
and associated parameter values that best fit the data. In 
that paper we focused on the problem of dealing with the 
multitude of potential models of die system by exploiting 
qualitative diagnosis techniques to generate a set of can- 
didate qualitative diagnoses, and we described two param- 
eter estimation techniques to deal with estimating the pa- 
rameters associated with the model, particularly when er- 
roneous behavior manifested itself some period of time 
after the initial occurrence of a fault (See (MBCGOO; 
MBCG99) for details.) We did not discuss the specific prob- 
lem of tracking multiple candidate models, nor did we dis- 
cuss how to compare them. 

In this paper, we formulate die hybrid monitoring and 
diagnosis task as a Bayesian model tracking and selection 
problem (e.g., (Mac91)). In particular, we wish to estimate 
the state (model) of the system at successive time instants, 
given a history of observations. The system diagnosis is do- 


scribed by the value of a specific subset of the state variables 
— namely those that designate whether components are nor- 
mal or abnormal, and what their associated parameter values 
are. We estimate state by tracking die posterior distribution 
of the state, given the observations. 

Probabilistic tracking of complex hybrid systems for diag- 
nosis purposes presents a number of interesting challenges. 
Kalman filtering techniques, traditionally used for tracking 
linear dynamical systems with Gaussian noise, assume a 
Gaussian density which is unimodal, making a Kalman fil- 
ter (Kal60) inadequate for simultaneously tracking alterna- 
tive candidate models. Multiple Kalman filters, one for each 
candidate model, can sometimes be used to track multiple 
candidate models of linear dynamical systems with Gaus- 
sian noise (e.g., (Fra90». More importantly, hybrid systems 
often have complex nonlinear, nonGaussian and potentially 
nondeterministic behavior. The nonlinearities come from 
both die mode switching (faulty or normal modes of behav- 
ior), and from the nonlinear dynamics within a mode. The 
latter has been addressed in some cases by using local lin- 
ear (Taylor series) approximations of the nonlinear contin- 
uous dynamics, such as is done with Extended Kalman Fil- 
ters (e.g., (BF88)) or Iterated Extended Kalman Filters (e.g., 
(Jaz70)). 

In this paper, following research on bootstrap filters, par- 
ticle filters and the condensation algorithm (e.g., (GSS93; 
IB98)), we use a factored sampling technique to sample and 
represent our multimodal posterior distribution of the state 
(models) given the observations. Such a technique enables 
us to track multiple models of nonlinear systems simulta- 
neously. Unfortunately, sampling techniques for probabilis- 
tic tracking focus on die most likely models within the dis- 
tribution, whereas most fault models have low probability, 
initially. To overcome this bias, we show how to integrate 
the qualitative diagnosis techniques described in (MBCGOO; 
MBCG99) into the temporal prior of our Bayesian formula- 
tion to focus sampling on models that are indicated by our 
qualitative candidate diagnoses. 

In die next section, we provide a brief description of 
NASA’s Sprint AERCam, which we have used as a motivat- 
ing example and which we will use to illustrate certain con- 
cepts in this paper. In die section that follows the description 
of the AERCam, we present a formal characterization of die 
class of hybrid systems we study and the diagnosis problem 
they present. Next, we describe our Bayesian formulation 
of die problem and the algorithm we use for computing and 
propagating posterior distributions. In the final section, we 
summarize, discuss our continuing research in this area, and 
reference some related work. 

The AERCam 

We are using NASA’s Sprint AERCam and a simulation of 
system dynamics and the controller written in Hybrid CC 
(HCC) (AG98) as a testbed for this work. To make this 
paper somewhat self-contained, we condense and repeat the 
description provided in (MBCG99). The AERCam is sim- 
pler than many of die complex systems we intend to diag- 
nose, but it serves well in illustrating the concepts developed 



velocities (u, v, w) are the components of the translation 
velocity, (p, q, r) are components of the angular velocity. 



Three views of the AERCam, showing the thrusters, 
and showing all the thrusters together in the cube 
circumscribing the AERCam. 

Figure 1: The AERCam axes and thrusters 


here, and has provided an excellent testbed for our prelimi- 
nary work. We describe the dynamic model of the AERCam 
system briefly, a more detailed description of the model and 
simulation appear in (AG98). 

The AERCam is a small spherical robotic camera unit, 
with 12 thrusters that allow both linear and rotational mo- 
tion (Fig. 1). For the purposes of this model, we assume the 
sphere is uniform, and the fuel that powers the movement is 
in the center of the sphere. The fuel depletes as the thrusters 
fire. 

The dynamics of the AERCam are described in the AER- 
Cam body frame of reference. The translation velocity of 
this frame with respect to the shuttle inertial frame of ref- 
erence is 0. However, its orientation is the same as the ori- 
entation of the AERCam, thus its orientation with respect 
to the shuttle reference frame changes as the AERCam ro- 
tates (i.e., it is not an inertial frame). The twelve thrusters 
are aligned so that there are four along each major axis in 
Ae AERCam body frame. For modeling purposes, we as- 
sume the positions of the thrusters are on die centers of die 
edges of a cube circumscribing the AERCam. Thus, for ex- 
ample, thrusters T \ . To ? T-j , T\ are parallel to the j*-axis and 
are used for translation along the r-axis or rotation around 



the .{/-axis. I.e., firing thrusters T\ and To results in transla- 
tion along the positive o*-axis, and firing thrusters T ( and T.\ 
results in a negative rotation around the y-axis. AERCam 
operations are simplified by limiting them to either transla- 
tion or rotation. Thrusters are either on or off, therefore, the 
control actions are discrete. In a normal mode of operation, 
only two thrusters are on at any time. 

AERCam dynamics 

A simplified model of the AERCam dynamics based on 
Newtonian laws is derived using an inertial frame of ref- 
erence fixed to the space shuttle. The AERCam position in 

this frame is defined as the triple (:r, y, r). Let V be the 
velocity in the AERCam body frame, with its vector compo- 
nents given by («, i\ w). The frame rotates with respect to 
the inertial reference frame with velocity o.’ — (p, ry, r), die 
angular velocity of the AERCam. The rotating body frame 
implies an additional Coriolis force acting upon the AER- 
Cam. We assume uniform rotational velocity since in the 
normal mode of operation, the AERCam does not translate 
and rotate at the same time (Am78, pg. 130). Similar equa- 
tions can be derived for the rotational dynamics (AG98). 

dim V)/dt =F — 2m(V xJ) Newton’s Law 

V dm /dt + md(y)/dt =F —2 m0 x V) 

The resultant equation for each coordinate: 

du/dt = Fj /m — 2 iqw — wr) — in/m) • dm/dt 

dv/dt = Fy/m — 2(rt* — pw) — (v/m) * dm/dt 

dw/dt = F ; /m — 2(jw — qu) — (w/m) * dm/dt. 

where the force F on each axis, is a function of the percent- 
age degradation of the thrusters that are exerting force in that 
direction as specified in Figure 1. Under normal operating 
conditions, the thrusters operate at 100 %. 

We use these models to predict the position of die AER- 
Cam at time / -J- 1, given the position at time L We add noise 
to each of die models above. In this case die noise is white 
Gaussian noise with a mean of zero and a standard deviation 
<7. As noted above, these models are implemented in HCC 
and are used to compute the likelihood described in the next 
section. 

Position Control Mode of the AERCam 

In the position control mode, die AERCam is directed to go 
to a specified position and point the camera in a particular di- 
rection. Assume the AERCam is at position A and directed 
to go to position B. In the first phase, the AERCam rotates 
to get one set of thrusters pointed towards B. These are then 
fired, and the AERCam cruises towards B. Upon reaching a 
position close to B, it fires thrusters to converge to B, and 
then rotates to point the camera in the desired direction. 

To facilitate the illustration of the diagnosis problem, we 
use a simple trapezoidal controller, which we explain in two 
dimensions. Suppose the task is to travel along the r-axis 
for some distance, then along the y-axis. Such manoeuvres 
are needed for navigating in the space shuttle. In order to do 


this, the AERCam fires its x thrusters for some time. Upon 
reaching the desired velocity, these are switched off. When 
the AERCam has reached a position close to the desired x 
position, the reverse thrusters are switched on, and the AER- 
Cam is brought to a halt — the velocity graph is a trapezium. 
The process is analogous for the y direction. 

Problem Formulation 

In this section we describe our formulation of the hybrid di- 
agnosis problem. Once again, the hybrid systems we ex- 
amine are discrete-time hybrid systems. Observations and 
state estimation are made at Tegular intervals 1.2,... , /, t + 
1. — Further, we assume that our systems contain no au- 
tonomous jumps. I.e., all nominal transitions between sys- 
tem modes are induced by a controller action, none are in- 
duced by the system state and mode (Bra95). Autonomous 
jumps are common in hybrid models where a mode with 
complex nonlinear behavior has been simplified by creating 
multiple modes of less complex behavior, with state-induced 
autonomous jumps connecting them. Building on tbe con- 
cepts in (MBCG00): 

Definition 1 (Hybrid System) A hybrid system is a 5-tuple 

(M, A\S,r,/>: 

• /t € M is the discrete state or mode of the system, where 
M is a finite collection of variables, /i/ is the system 
mode at time / . 

• x € X C R n is the continuous state vector of the system. 
Xf is die continuous state at time / . 

• (T € E, is the discrete input, where E is a finite collection 
of actions. I.e., the controller actions that transition the 
system between modes. 

• v € V C if" is the continuous input. 

• / is the system dynamics function that maps the mode, the 
continuous state, and the input into the mode and contin- 
uous state at the next discrete time point. (/i /+I , .r /+ ! ) = 
/(/«/, x-i , 07 , V/ . w f ), where wt € R m is zero-mean white 
noise of known pdf, and / : M x .Y x E x R v x R m -* 
M x A'. / is often expressed as a collection of func- 
tions, e.g., functions that describe the continuous behav- 
ior within a specific mode, and a function that describes 
the discrete transitions between modes, based on discrete 
input. 

• d)s € R v is the observation vector of the system. obs { is 
the observation vector at time I. ot>s ( is related to the con- 
tinuous state vector x t by the function oft*, — h(x t . u f ) 
where v t € R r is zero-mean white noise of known pdf, 
and A: A' x R r -4 W'. 

Definition 2 (System State) The state of a hybrid system at 
time lj {pi , Xf ) comprises the discrete mode of die system 
and the continuous state at L 

To define the hybrid diagnosis problem, we augment Defini- 
tion 1 as follows. 

Definition 3 (Diagnosable Hybrid System) A diag- 
nosable hybrid system , (M. X. E, V. /. COMPS) is a 
hybrid system comprised of rn potentially malfunctioning 
components COM PS = ,c m ) where 



• For each ft € M, /< includes a designation of whether 
each a € COMPS is operating normally, or abnormally, 
i.e., 

• For each ft, continuous state vector r includes a set of 
distinguished parameters 0 associated with that mode. 

• We assume that transitions to fault modes are achieved by 
exogenous actions. Hence, E = E r U E f , where 

- E r is a finite set of controller actions, and 

- E r is a finite set of exogenous actions. 

We introduce the following additional notation, 

• O, designates the observation history, the sequence of 
time-indexed observations. Ot designates the observation 
histoiy to time L 

• ftp denotes a faulty mode, i.e., a mode for which at least 
one Ci € COMPS is ab(ci) in ftp . Op denotes the pa- 
rameters associated with ftp. 

In the case of the AERCam example, the potentially mal- 
functioning components are the 12 thrusters, and a mode 
ft includes die behavior mode (e.g., translate-x, translate- 
y, rotate-x, etc.) and [~]ab(Ti),i — 1 , ... .12, for each 
thruster. The continuous state vector includes the -r, t/, r 
position of the AERCam, velocity and acceleration. The pa- 
rameter values, 0 associated with each ft are the percentage 
degradation of each of the thrusters. As we will see later 
on, we make a Maikov assumption with respect to comput- 
ing the temporal dynamics of our system. Hence all relevant 
state must be included explicitly in the state variables. 

Definition 4 (Model) A model of a diagnosable hybrid sys- 
tems is a time-indexed mode sequence and associated pa- 
rameter values ([ft i , . . . , ft m ]* [<9| . . . . , #,„]). The model to 
time l is denoted (//. 0) and the model at time / is denoted 
s t = (ftj.Oi). The model is a distinguished subset of the 
entire system state. 

In this paper we make several simplifying assumptions re- 
garding our diagnosis task. In particular, we make a single- 
time fault assumption. We assume that our systems do not 
experience multiple sequential faults. Further, we assume 
that faults are abrupt, resulting in partial or full degradation 
of component behavior. We cast the hybrid diagnosis task 
as the problem of finding the most likely model for the ob- 
servation history, P(s/ j O), i.e, the mode and parameter 
values (fit , Of ) that best fit the observations over time. To do 
this, we appeal to a Bayesian formulation of the problem. 

Bayesian Formulation 

To monitor and diagnose a hybrid system, we must compute 
the posterior probability distribution over models at time /, 
given die observation histoiy. Recall, using Bayes’ rule that 
the posterior is proportional to the likelihood times the prior. 
I.e., 

p(mode\ | observations) oc p ( observations | model) (model). 

Our objective is to find the posterior probability distribu- 
tion over models at time /, s t , given the observation history 
up to time t, Ot . I.e, we wish to compute p(s f | Ot ) . 


To compute the temporal dynamics of our system, we 
make a Maikov assumption, i.e., 

p(s, | a/-!,... .* 0 ) =p(*/ | si- i) 

Further, we assume that at each time point, there is a small 
probability of an exogenous action, leading to a transition 
to a failure mode. Finally, we assume that given the current 
model .s/, the current observations obs t and previous obser- 
vation histoiy Ot-\ are independent. 

Hence, in order to track our hybrid system, we can com- 
pute the posterior distribution of the model at time l given 
the observation histoiy which, according to Bayes’ rule 
and our assumptions above, is proportional to the likeli- 
hood of the observation at time / given the model at time 
t (p(ol)Si | .$/)) and the temporal prior, the prediction 
of the current model, given the observation history up to 
/ - l,(p(*/ | Ot- i). I.e., 

p(s, | O t ) “ kp(d)S t | s t ) p(st | Ot- 1 ), 

where k ensures that die distribution integrates to one. 

The likelihood of the observations given the state is easily 
evaluated for the AERCam following the model described in 
die previous section. The temporal prior, i.e., die probability 
of the current model given the observation history to / - 1 
depends on the posterior over models at the previous time 
point, p(st - 1 | <3/_i) and the temporal dynamics, p(Sf \ 
I.e., 

p(*i | 0,-i ) = f p(si | .s,_i)p(.s,_i | 0,-i Ms, -I 
J *i - 1 

The temporal prior expresses the probability of a partic- 
ular model given the observation history up to that point. 
In the case of a fault diagnosis, the likelihood of a fault 
model will initially be very low. If we are tracking using a 
finite number of parallel filters, or using a factored sampling 
method as suggested in the next section, this may mean that 
we will initially not track these fault models, or alternately 
that we track many low probability models which is com- 
putationally expensive. In order to focus the temporal prior 
more quickly and accurately on the appropriate diagnostic 
models, we make use of qualitative diagnosis techniques. 

In (MBCGOO; MBCG99), we proposed to use qualitative 
diagnosis techniques to generate qualitative candidate diag- 
noses - candidate mode and parameter values that were con- 
sistent with observations O in some window of time. 

Definition 5 (D-tuple (MBCGOO)) A D-tuple is a 4-tuple 
(C. ftp Jr -Op), where ftp is a fault mode, tp is the time 
file fault mode commenced. Op is the parameter values as- 
sociated with tiie fault mode behavior, and C is the set of 
failed (abnormal) components in ftp. 

Definition 6 (Candidate Qualitative Diagnosis (MBCGOO)) 

Given a diagnosable hybrid system with model ( ft, 0), in- 
put history J 1 , and observation histoiy, O, D-tuple 
(C. ft pjp ? Op) is a candidate qualitative diagnosis iff there 
exists a range of parameter values Op = [Oi, O n ], and time 
range Ip = [//. t u ] such that the occurrence of fault mode 
ftp with parameter values Op in time range t p is consistent 
with O, X and (/7. §). 



We do not repeat the diagnosis algorithms here, but re- 
fer the reader to (MBCGOO; MBCG99) for details. These 
generated diagnoses are used to propose a set of different 
models to be tracked by the system. The candidate models 
are generated by exploiting previous work on qualitative di- 
agnosis of continuous systems (e.g., (MB99)), adapting the 
authors' causal propagation algorithms to deal with the dis- 
crete state variables and mode transitions of the hybrid sys- 
tems. To incorporate this so-called oracle into our Bayesian 
formulation, we use it to bias or focus the temporal prior. 
This will in turn more heavily weight die posterior for the 
corresponding fault models, s/ . In the case of particle filter- 
ing, the technique we propose in the next section to compute 
die posterior, this focusing of the temporal prior will help 
die algorithm sample from the appropriate part of the dis- 
tribution. To incorporate this qualitative diagnosis “oracle” 
we may alter our view of the posterior we are computing as 
follows. 

p(s f | O, . oracle) « p(ofts f | s t , oracle) 

p(.sv | (9/ _i. oracle) 
a p{obfit | */) p{Xf | Of~t , oracle) 

where p(s t \ Oi -i, oracle) is equal to p(s f | (9,_i) above, 
when the observations are consistent with the current model, 
and otherwise p(x f | <9/_ j , oracle) is simply the normalized 
probability of the faulty models, given the observations. To 
ensure the speed of the oracle, and because of the lack of re- 
liable numbers for such calculations, the probabilities gener- 
ated by the oracle are normalized prior probabilities of dif- 
ferent faults given the observations, as defined by die system 
builder. 

Once the posterior is computed, different models can be 
compared by estimating die expected value of different mod- 
els, normalizing and comparing. For example, we may sum 
the likelihoods for all samples having like [->]a6(c f ) desig- 
nations, and compare these to determine which components 
are likely malfunctioning. 

Computing the Posterior 

In die previous section we presented the problem of tracking 
and diagnosing hybrid systems using a Bayesian formula- 
tion. As noted in die introduction, there are many algorithms 
for probabilistic tracking of dynamical systems, though most 
are not tailored to simultaneously tracking multiple candi- 
date models nor to dealing with nonlinear dynamics. Our 
posterior distribution p(s t | Q \ ) will be a multi-dimensional, 
multi-modal distribution, reflecting the multiple competing 
diagnostic models. There is no closed-form (parametric) 
representation for this distribution, as there is, for exam- 
ple, for a unimodal Gaussian. Consequently, to compute this 
posterior, we appeal to factored sampling techniques to pro- 
vide an approximation of the distribution, and project this 
distribution forward through time according to its dynamics, 
using the Condensation algorithm (IB98), derivative of the 
bootstrap algorithm (GSS93) and commonly referred to as a 
particle filter. 

1 Previously referred to as the action history. 


More specifically, die posterior distribution p(s t \ (9/), 
is represented as a set of N weighted samples {s ( 1 >, ... , 
with associated weights {7r (l) ,... ? 7r (>f) }. Intu- 
itively, the larger the N, the better the approximation, but 
the more costly the computation. Hence we would like to 
sample die distribution as sparsely as possible, while maxi- 
mizing our coverage of our distribution, and thus weighting 
samples more heavily in those parts of the distribution that 
have greater volume. 

At each time step, the basic algorithm comprises three 
steps: select, predict, and update. 

Select: We start with the posterior from the previous time 
step, p(*/_ i | Ot-\), represented as the factored sample 
(sj^ , , it}!!-, ), i = 1, . . . , N. Sample N times with replace- 
ment with probability 7r}^, , the sample {sj^ ( ), producing 

the samples {s'} 0 }. Note that samples with high weights 
may be chosen multiple times. 

Predict: For each new sample s'}*’, propagate the sample 
forward according to the dynamics of the system to pro- 
duce new samples {s} l) ). In the case of our AERCam, these 
are the dynamics described in the previous section, together 
with zero-mean Gaussian white noise. This new set of sam- 
ples approximates a fair random sample for die effective 
prior p(s/ | (9/_ i ). What remains to compute is the weights. 
Update: Compute the weights, tt^ = p{ofjx ( \ s f — sj #) ). 
From the observations ol>s t , evaluate the likelihood of each 
sample, and normalize the likelihoods of the samples so they 
sum to 1. 1.e., 


I si 0 ) 

£»=t I s) n) ) 


The above algorithm does not reflect our qualitative diag- 
nosis oracle. In order to suitably focus the temporal prior, we 
use a linear combination of the samples from the computed 
temporal prior, and samples from die oracle. This technique 
was inspired by (BF99), and could also be achieved using 
importance sampling. 

The sample approximation to the distribution, p(s t \ Oi) 
can be used to compute the expected value for some moment 
/ of the density, for example a mean of some state variable, 
i.e., 


7—1 

In this way, we can compare the sum of the likelihoods for 
each distinct model. 

Summary and Related Work 

In this paper we expanded the hybrid diagnosis framework 
described in (MBCG99; MBCGOO) to present a mathemati- 
cal formulation and computational techniques for generating 
diagnoses of hybrid systems in terms of Bayesian tracking 
and model comparison. We characterized the evaluation of 
our models (system mode and associated parameter values) 



as the computation of the posterior distribution of models 
given a history of observations. Exploiting a Markov as- 
sumption, we showed that this could be computed in terms 
of die likelihood of the observations at time /, given the 
model at time /, times a prior. Exploiting the work de- 
scribed in (MBCG99; MBCGOO) for generating qualitative 
diagnoses of hybrid systems, we treated our qualitative mon- 
itoring and diagnosis system as an oracle. If the observations 
were consistent with die current model, then the qualitative 
monitoring and diagnosis system had no effect on the com- 
putation of the posterior. However, if the observations were 
inconsistent then the oracle would generate a set of candi- 
date diagnoses that would be used to adjust the prior to focus 
the likelihood computation on that part of the model space 
that was indicated by the qualitative monitoring and diagno- 
sis engine. 

Since hybrid systems are generally nonlinear, and hence 
the distribution of the posterior multimodal and non- 
Gaussian, we represented the posterior distribution as dis- 
crete samples and exploited factored sampling techniques, 
used in particle filtering and in the Condensation algorithm, 
to propagate conditional probability densities over time. 

We are still in die early stages of experimenting with these 
techniques, but preliminary results look promising. Con- 
densation has proven effective for some near realtime visual 
tracking tasks (e.g., (IB98)), but we anticipate that more 
complex hybrid systems with large state spaces and par- 
tial observability will require further computation and larger 
amounts of memory that will compromise realtime compu- 
tation, just as they do, for example, with POMDPs. Such 
systems will require new variants of many of the techniques 
we currently employ in model-based diagnosis including ex- 
ploiting problem decomposition, compact representations of 
state spaces, abstractions of problems, and approximation of 
inference. In summary, Bayesian tracking and model com- 
parison and factored sampling techniques for dynamical sys- 
tems provide a sound mathematical formalism and promis- 
ing tools for monitoring and diagnosing complex dynamical 
systems. 

The problem of monitoring and diagnosing hybrid sys- 
tems has received little attention to date, although there is 
much related work. Within the AI community, there has 
been a great deal of research on diagnosing static systems 
(e.g., (HCD92)), while much less on diagnosing discrete dy- 
namical systems (e.g., (CT94; McI98; WN96; BLPZ99)), 
qualitative diagnosis of continuous systems (e.g.,. (MB99)), 
and tracking (e.g., (RK99)). Most recently, (LPKBOO), have 
developed related techniques for monitoring and diagnosing 
Conditional Linear Gaussian hybrid systems using a Dy- 
namic Bayes Nets to compactly represent the conditional 
probability distribution, and proposing algorithms for hy- 
pothesis reduction and smoothing. Within the FDI commu- 
nity, the largest proportion of research has focused on diag- 
nosing continuous systems (e.g., (Ger98; Fra90)). These ap- 
proaches have often used observer schemes and/or Kalman 
filters to track continuous system behavior. Diagnosis of 
discrete-event systems has also been studied within the FDI 
community (e.g, (SSLST96; Lun99)). Nevertheless, our 
work and the concurrent work of (LPKBOO) has been Ihe 


first to propose a Bayesian tracking approach to diagnosing 
hybrid systems. Our use of factored sampling techniques 
and particle filtering drawn from ft le statistics and computer 
vision c ommuni ties, presents a significant contribution to a 
challenging problem. 
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Abstract. This paper reports on an on-going project to investigate techniques to 
diagnose complex dynamical systems that are modeled as hybrid systems. In par- 
ticular, we examine continuous systems with embedded supervisory controllers 
that experience abrupt, partial or full failure of component devices. We cast the 
diagnosis problem as a model selection problem. To reduce die space of potential 
models under consideration, we exploit techniques from qualitative reasoning to 
conjecture an initial set of qualitative candidate diagnoses, which induce a smaller 
set of models. We refine these diagnoses using parameter estimation and model 
fitting techniques. As a motivating case study, we have examined the problem of 
diagnosing NASA’s Sprint AERCam, a small spherical robotic camera unit with 
12 thrusters that enable both linear and rotational motion. 

1 Introduction 

The objective of our project has been to investigate how to diagnose hybrid systems 
- complex dynamical systems whose behavior is modeled as a hybrid system. Hybrid 
models comprise both discrete and continuous behavior. They are typically represented 
as a sequence of piecewise continuous behaviors interleaved with discrete transitions 
(e.g., [7]). Each period of continuous behavior represents a so-called mode of the sys- 
tem. For example, in the case of NASA’s Sprint AERCam, modes might include trans- 
late JC-axis, rotate JC-axis, translate. Y-axis, etc. [ I], In the case of an Airbus fly-by-wire 
system, modes might include take-off, landing, climbing, and cruise. Mode transitions 
generally result in changes to the set of equations governing the continuous behavior of 
the system, as well as to the state vector that initializes that behavior in the new mode. 
Discrete transitions that dictate mode switching are modeled by finite state automata, 
temporal logics, switching functions, or some other transition system, while continuous 
behavior within a mode is modeled by, e.g., ordinary differential equations (ODEs) or 
differential and algebraic equations (DAEs). 

The problem we address in this paper is how to diagnose such hybrid systems. For 
the purposes of this paper, we consider the class of hybrid systems that are continuous 
systems with an embedded supervisory controller, but whose hybrid models contain no 
autonomous jumps. I.e., all nominal transitions between system modes are induced by 
a controller action, none are induced by the system state and model [7]. The class of 
systems we consider can be modeled as a composition of a set of component subsys- 
tems, each of which is itself a hybrid system. We assume that the system operation is 
being tracked by a monitoring and observer system (e.g., [ 19]) that ensures that the sys- 
tem behavior predicted by the model does not deviate significantly from file observed 


behavior in normal system operation. When observations occur outside this range, the 
behavior is deemed to be aberrant and diagnosis is initiated In this paper, we consider 
faults whose onset is abrupt, and which result in partial or complete degradation of 
component behavior. The general problem we wish to address can be stated as follows: 
Ci ve/i a hybrid model of system behavior, a history of executed controller actions, a his- 
tory of observations, including observations of aberrant behavior relative to the model, 
isolate the fault that is the cause for the aberrant behavior. Diagnosis is done online 
in conjunction with the continued operation of the system. Hence, we divide our diag- 
nosis task into two stages, initial conjecturing of candidate diagnosis and subsequent 
refinement and tracking to select the most likely diagnoses. 

In this paper we conceive the diagnosis problem as a model selection problem. The 
task is to find a mathematical model and associated parameter values that best fit the sys- 
tem data. These models dictate the components of the system that have malfunctioned, 
their mode of failure, die estimated time of failure and any additional parameters that 
further characterize the failure. To address this diagnosis problem, we propose to ex- 
ploit AI techniques for qualitative diagnosis of continuous systems to generate an initial 
set of qualitative candidate diagnoses and associated models, thus drastically reducing 
the number of potential models for our system. This is followed by parameter estima- 
tion and model fitting techniques to select the most likely mode and system parameters 
for candidate models of system behavior, given both past and subsequent observations 
of system behavior and controller actions. The main contributions of the paper are: 1) 
formulation of the hybrid diagnosis problem; 2) the exploitation of techniques for qual- 
itative diagnosis of continuous systems to reduce the diagnosis search space; and 3) the 
use of parameter estimation and data fitting techniques for evaluation and comparison 
of candidate diagnoses. 

In Section 2 we provide a brief description of NASA’s Sprint AERCam, which we 
have used as a motivating example and which we will use to illustrate certain concepts 
in this paper. In Section 3 we present a formal characterization of the class of hybrid 
systems we study and the diagnosis problem they present. In Section 4 we describe our 
approach to hybrid diagnosis and the algorithms we use to achieve hybrid diagnosis. 
Hie generation of initial candidate qualitative diagnoses is described in Section 4.1, 
and the subsequent quantitative fitting and tracking of candidate diagnoses and their 
models is described in Section 4.2. In the final two sections, we briefly discuss related 
work and summarize our contributions. 

2 Motivating Example: The AERCam 

We are using NASA’s Sprint AERCam and a simulation of system dynamics and die 
controller written in Hybrid CC (HCC) as a testbed for this work. We describe the 
dynamic model of the AERCam system briefly, a more detailed description of the model 
and simulation appear in [ 1] . 

The AERCam is a small spherical robotic camera unit, with 12 thrusters that allow 
both linear and rotational motion (Fig. 1). For the purposes of this model, we assume 
the sphere is uniform, and the fuel that powers the movement is in die center of the 
sphere. The fuel depletes as die thrusters fire. 



The Body frame of reference 
and the directions of velocities 
(u,v,w) are the components of 
the translation velocity, while 
(p,q,r) are components of the 
angular velocity. 



Three views of the AERCam, showing the thrusters, 
and showing all the thrusters together in the cube 
circumscribing the AERCam. 



Fig. 1. The AERCam axes and thrusters 


The dynamics of the AERCam are described in the AERCam body frame of refer- 
ence. The translation velocity of this frame with respect to the shuttle inertial frame of 
reference is 0. However, its orientation is the same as the orientation of the AERCam, 
thus its orientation with respect to the shuttle reference frame changes as the AERCam 
rotates (i.e., it is not an inertial frame). The twelve thrusters are aligned so that there 
are four along each major axis in the AERCam body frame. For modeling purposes, 
we assume the positions of the thrusters are on the centers of the edges of a cube cir- 
cumscribing the AERCam. Thus, for example, thrusters Ti,T*,Tz,Ti are parallel to 
the j*-axis and are used for translation along the :r-axis or rotation around the //-axis. 
I.e., firing thrusters T\ and X ' 2 results in translation along the positive :r-axis, and firing 
thrusters T\ and T\ results in a negative rotation around the y- axis. AERCam operations 
are simplified by limiting them to either translation or rotation. Thrusters are either on 
ot off, therefore, the control actions are discrete. In a normal mode of operation, only 
two thrusters are on at any time. 


2.1 AERCam dynamics 

A simplified model of file AERCam dynamics based on Newtonian laws is derived us- 
ing an inertial frame of reference fixed to the space shuttle. The AERCam position in 

this frame is defined as the triple (a\ y. :). Let V be the velocity in the AERCam body 
frame, with its vector components given by («, v, w). The frame rotates with respect 
to the inertial reference frame with velocity u> = ( p , q \ r), the angular velocity of the 
AERCam. The rotating body frame implies an additional Coriolis force acting upon the 
AERCam. We assume uniform rotational velocity since in the normal mode of opera- 



tion, the AERCam does not translate and rotate at the same time [2, pg. 130]. Similar 
equations can be derived for the rotational dynamics [1]. 

dim V)/dt — I? — 2m(V x j) Newton’s Law 
V dm/ dt + md(y)/dt =T -2m(5 x V) 

The resultant equation for each coordinate: 

du/dt = F x /m — 2{qw — vr) — (w/m) * dm/dt 
dvjdt = F,,/m — 2 (ru — ptr) — iv/m) * dm/dt 
dw/dt = Ft/m — 2 ijm — qu) — (w/m) * dm/dt 


2.2 Position Control Mode of the AERCam 

In die position control mode, die AERCam is directed to go to a specified position and 
point the camera in a particular direction. Assume the AERCam is at position A and 
directed to go to position B. In the first phase, the AERCam rotates to get one set of 
thrusters pointed towards B. These are then fired, and the AERCam cruises towards B. 
Upon reaching a position close to B, it fires thrusters to converge to B, and then rotates 
to point the camera in die desired direction. 

To facilitate the illustration of the diagnosis problem, we use a simple trapezoidal 
controller, which we explain in two dimensions. Suppose the task is to travel along 
the r-axis for some distance, then along the </- axis. Such manoeuvres are needed for 
navigating in the space shuttle. In order to do this, the AERCam fires its ;r thrusters 
for some time. Upon reaching the desired velocity, these are switched off. When the 
AERCam has reached a position close to the desired ;r position, the reverse thrusters are 
switched on, and die AERCam is brought to a halt — the velocity graph is a trapezium. 
The process is analogous for the y direction. 


3 Problem Formulation 

In this section we provide our formulation of the hybrid diagnosis problem. 

Definition 1 (Hybrid System). A hybrid system is a 5-tuple {M, X. T, X, <6), where 

- M y finite set of system modes (yu , . . . , /<*.). 

- A' C R n , continuous state variables. :r(l) is the continuous behavior at time / . 

- Ty finite set of functions {/ /t1 , . . . , f fth }, and associated parameter values 0 such 
that for each mode, //,, / M( (/.^,:r(/)) : R x R x X -> X defines the continuous 
behavior of the system in fu? 

- X, finite set of actions (<r \ , . . . ? <t/), which transition the system between modes. 

- </>, transition function which maps an action, mode and system state vector into a 
new mode and initial state vector, i.e., <j> : X x M x X -> M x A'. 

To define the hybrid diagnosis problem, we augment Definition 1 as follows. 

1 Parameter value ranges may be associated with Q. 



Definition 2 (Diagnosable Hybrid System). A diagnosable hybrid system, 
{M r X , f, X. 0 , COMPS) is a hybrid system comprised of m potentially malfunc- 
tioning components COMPS = [r\ ..... r m ) where 

- For each fi € M y includes a designation of whether each r, 6 COMPS is 
operating normally, or abnormally, i. e., ( ) . 

- We assume that transitions to fault modes are achieved by exogenous actions. 
Hence, X = X r U X e , where 

• X f is a finite set of controller actions, and 

• X e is a finite set of exogenous actions. 

- A, the controller action history, the sequence of time-indexed controller actions 
performed. 

- Xou* Q X , continuous state variables that are observable. :r o/ , 5 (/) is the observa- 
tions at time /. 

- O, the observation history, file sequence of time-indexed observations. 

For notational convenience, ftp denotes a faulty mode, i.e., a mode for which at least 
one Ci € COMPS is ab(ci) in ftp. 9p denotes the parameters associated with f flF . 

In file case of the AERCam example, the potentially malfunctioning components are 
the 12 thrusters, and a mode ji includes the behavior mode (e.g., translate-x, translate- 
y, rotate-x, etc.) and (~> )ab(Ti). 1 = 1, . . . . 12, for each thruster. The continuous state 
vector includes the :r, y, z position of the AERCam, velocity and acceleration. The 
parameter values, 6 associated with each f it are the percentage degradation of each of 
the thrusters. 

Definition 3 (Model)* A model, M od of a diagnosable hybrid systems is a time-indexed 
mode sequence and associated parameter values ([/t j ..... //„,]. [9\ ..... 9 tn ]) 

Notice that each model of the system, (/*, 0) induces a corresponding time-indexed 
piecewise continuous sequence of functions [/ /t) , . . . . ] dictating system behavior. 

In this paper we make several simplifying assumptions regarding our diagnosis task. 
In particular, we make a single-time fault assumption. We assume that our systems do 
not experience multiple sequential faults. Further, we assume that faults are abrupt, 
resulting in partial or full degradation of component behavior. We cast the hybrid diag- 
nosis task as the problem of finding the most likely model for the observation history, 
P(Mod | <f>). I.e, the sequence of modes and parameter values (/*. 9) that best fit the 
observations overtime. Under normal operation, the model of the system Mod nornu j is 
fully dictated by the sequence of controller actions -4 and the nominal parameter values, 
9 . Once again, we assume that the system operation is being tracked by a monitoring and 
observer system (e.g., [ 19]) that ensures that the system behavior predicted by the model 
does not deviate significantly from the observed behavior in normal system operation. 
When observations occur outside this range, the behavior is deemed to be aberrant and 
diagnosis is initiated. Given a diagnosable hybrid system (M, X . T, X , <?, COMPS), 
a controller action history, A and a history of observations, O which includes observa- 
tions of aberrant behavior, file hybrid diagnosis task is to determine what components 
are faulty, what fault mode caused die aberrant behavior, when it occurred, and what the 
values of the parameters associated with the fault mode are. In the AERCam system, a 
diagnosis might be that thruster Tj experienced a blockage fault of 50%, at time / 



Once Mod nnrmti i has been rejected, we must find a new most likely model from 
among the potentially exponential (in C OM PS) number of mode sequences, occurring 
within a large but bounded time range. We propose to exploit previous research on 
temporal causal graphs for qualitative diagnosis of continuous systems [ 1 8], to compute 
a set of candidate qualitative diagnoses that are consistent with our system, in order to 
identify a preliminary subset of candidate models, whose likelihood can be estimated. 

Definition 4 (D-tuple), A D-tuple is a 4-tuple (C,fipJp,0p) 9 where ftp is a fault 
mode, Ip is the time the fault mode commenced. Or is the parameter values associated 
with the fault mode behavior, and C is the set of failed (abnormal) components in ft p. 

Definition 5 (Candidate Qualitative Diagnosis). Given a diagnosable hybrid system 
with model Mod = (ft, 0) an action history A, and a history of observations, O which 
includes observations of aberrant behavior, D-tuple (C. ftpAr , Op) is a candidate qual- 
itative diagnosis iff there exists a range of parameter values Op = [9pO n ], and time 
range / r = [// , /„] such that the occurrence of fault mode ft p with parameter values Op 
in time range Ip is consistent with O, .4 and Mod. 

Hence, a candidate qualitative diagnosis stipulates a fault mode, including one or 
more faulty components. It also stipulates a lower and upper bound, [//, /•„], on the time 
the fault mode occurred. This range generally corresponds to the start times of the con- 
troller induced modes preceding and following the fault, or up to the point tire fault was 
detected. This candidate diagnosis induces an associated candidate model , Mode — 
(D* i IH, /* f , ..... /<!„], ,0i, Or, 0 ' i+ ...... O corresponding to M od 

with the fault mode ft p and Op inserted at i p. Every subsequent mode, , . . . , ft m , 
has nb(Ci), Ci € C enforced, and every subsequent set of parameters has the param- 
eters associated with faulty components C enforced. Computing candidate qualitative 
diagnoses is discussed in Section 4.1 . 

Since each candidate qualitative diagnosis only conjectured ranges for the time of 
the fault mode, / r and parameter values associated with the fault mode, Op, the asso- 
ciated candidate models are underconstrained. In Section 4.2, we discuss methods for 
estimating unique values for t p and Op and for estimating a posterior probability for 
each of tire candidate models. Mode , given O. 

Definition 6 (Candidate Diagnosis). Given a diagnosable hybrid system, a history of 
controller actions A, and a history of observations O , D-tuple (C, ftp, Ip, Op) with 
associated model Mode is a candidate diagnosis for the hybrid system, iff P(M ode I 
O) > o , for defined threshold value a € [0, 1]. 


4 Diagnosing Hybrid Systems 

In this section we discuss one method for computing hybrid diagnoses. In Section 4.1 
we discuss a technique for generating candidate qualitative diagnoses, and their associ- 
ated candidate models. In Section 4.2 we discuss techniques for model fitting and for 
model (and hence diagnosis) comparison. In particular we discuss techniques for esti- 
mating the parameters of the candidate models, and the likelihood of file models, and for 


continued monitoring and refinement of die candidate models as the system continues 
to operate and observations continue to be made. 

We illustrate these techniques with die following simple AERCam example. Con- 
sider the scenario depicted in Fig. 2. In the first accelerate phase, the AERCam is being 
powered by thrusters T 1 and T2. Assume that at some point in this phase, a sudden leak 
in the T2 thruster causes an abrupt change in its output. As a consequence, the AER- 
Cam starts veering to the right of die desired trajectory, as illustrated by the left-most 
dotted lines in Fig. 2. (The other dotted lines represent other potential candidate diag- 
noses consistent with the point of detection of the failure.) Soon after this occurs, the 
supervisory controller commands the AERCam to turn off Thrusters T 1 and T2 with 
the objective of getting the AERCam to cruise in a straight line. In the faulty situation, 
the AERCam has some residual angular velocity about die z-axis, so it continues to 
rotate in the cruise mode. Then the controller turns on thrusters T3 and T 4, to decel- 
erate the AERCam with the objective of bringing it to a halt. Again, this objective is 
not entirely achieved in die die faulty situation. Next, thrusters To and T6 are switched 
on, to move the AERCam in the y direction. However, since the AERCam is not in the 
desired orientation after the failure, the position error due to faulty thruster T2 accumu- 
lates causing a greater and greater deviation from the desired trajectory of the system. 
The position of the AERCam is being continuously sensed, filtered for noise and mon- 
itored. At some point within die y translation the trajectory exceeds the error bound, 
i.e., P(Mod norma i < a) and is flagged by the monitoring system as aberrant relative 
to Mod n0 r m ab At this point, the diagnosis task begins. 
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Fig. 2. Possible fruit trajectories of AERCam ( simplified for illustration purposes). 


4*1 Qualitative Candidate Generation 

Given the current system model Mod — (/*, 0) (commonly Mod normi j), a history of 
controller actions „4, and a history of observations including one or more observa- 


tions of aberrant behavior, we wish to generate a set of candidate qualitative diagnoses 
{C. ftp. Ip, Op), and associated candidate models as described in Definition 5. To do 
so, we extend techniques for generating qualitative diagnoses of continuous dynamic 
systems to deal with hybrid systems with multiple modes. The model and propagation 
mechanism, as applied to continuous systems diagnosis, is described in [18]. 

In the case of our AERCam example, the action history ^4 is [(on(Tl), on(T2)), 
(ofiRXl), offlT2)) ? (on(X3), on(T4)), (off(T3) ? off(T4) ? on(T3), on(TG)), (off(Xo), 
off[TG))]; the model, Mod nornu u is the time-indexed sequence [(accelerate jr, ^ab(T 1- 
T12), 0 ), (cruise^r. ^ab(Tl~TV2), 0),(deceleralejr 7 -iah(Tl-T12), 0), (accelerate -if, 
-^ah(Tl - T12), 0) y {cruise4i. -iab(Tl - T12), 0)], where 0 is a vector of length 12 all 
of whose entries are 0 (percent degradation in thrusters). 

To generate candidate qualitative diagnoses we construct an abstract model of the 
dynamic system behavior, Mod nwmfi i as a temporal causal graph. A part of the tem- 
poral causal graph for the AERCam dynamics is shown in Fig. 3. The graph expresses 
directed cause-effect relations between component parameters and die system state vari- 
ables. Links between variables are labeled as: (i) + 1, implying direct proportionality, 

(ii) - 1, implying inverse proportionality, and (iii) /, implying an integrating relation. 

An integrating relation introduces a temporal delay in that a change on the cause side of 
the relation affects the derivative of the variable on the effect side. This adds temporal 
characteristics to the relations between variables. Some edges are labeled by variables, 
implying the sign of die variable in the particular situation defines die nature of the rela- 
tionship. The candidate generation algorithm is invoked for every initial instance of an 



Fig. 3. A subset of die temporal causal graph showing the relations between Thrusters Tl — T8 
and the x and y positions of the AERCam. 


aberrant observation. The aberrant observation plus the controller action history .4 are 
input to a backward propagation algori thm that operates on the temporal causal graph. 



The algorithm operates backwards from the last mode in the mode sequence of M od: 

Step 1 For the current mode, extract the corresponding temporal causal graph model, 
and apply die Identify Possible Faults algorithm. Details of this algorithm are presented 
in [18], but the key aspect of this algorithm is to propagate the aberrant observation ex- 
pressed as a ± value, backward depth-first through the graph. For example, given that 
the ^-position of the AERCam has deviated - (i.e., below normal), backward prop- 
agation implies d(y)/dt is and so on, till we get Tf and Tf y implying thrusters 
To and TG are possibly faulty with decreased thrust performance. Propagation along a 
path can terminate if conflicting assignments are made to a node. The goal is to system- 
atically propagate observed discrepancies backward to identify all possible candidate 
hypotheses that are consistent with die observations. In our example, the component 
parameters, COMPS = {IT, . . . , T12) form the space of candidate faults. 

Step 2 Repeat Step 1 for every mode in the mode sequence, to i . The system model 
needs to be substituted as the algorithm traverses the mode sequence backwards. There- 
fore, back propagation will be performed on a different temporal causal graph for each 
mode in the controller history 2 . 

The output of this step is a set of qualitative diagnoses (C, f irJ p. #r}, each with 
an associated candidate model, as described in Section 3. Returning to our AERCam 
example, three qualitative candidate diagnoses are generated. The first candidate diag- 
nosis is that T2 failed in the x acceleration phase. The time of the fault mode transition 
is {/ 1 , h\, and the parameters associated with the failure - the percentage degradation 
of the component is in the range [0, 100], So the first candidate qualitative diagnosis 
is { T2 , (accelerate Jt. ab(T2), ^ab(T\.TZ - TT2), Op). [/ 1 , /•>], [0 ? 100]). The candi- 
date model simply has (accelerate ab(T2), ~^ab(Tl). ->ab(T3-T 12)) insertedafter 
the mode (accelerate sv. -» ab(Tl - T12)), and ab(T2) enforced in every subsequent 
mode. The second candidate qualitative diagnosis is that T4 failed in the deceleration 
phase of x translation, i.e., (T 4, (decelerate jr, ab{T 4) , ->ab(Tl - T3, To - T12) , 0 p). 
p 3 , /. t], [0, 100]). The third candidate is that TG foiled during y acceleration, i.e., (TG, 
( aecelcrate-y , ab(T6 ), ~*ab(Tl - Tb,T7 - T12),0p). [f.i, Ip], [0, 100]), where tp is 
the time of detection of the aberrant behavior. In each case Op is a vector of length 12 
with every entry equal to 0 (percentage degradation), except the entries corresponding 
to the faulty thrusters, C which will have the range [0, 100]. 


4.2 Model Fitting and Comparison 

Given the candidate qualitative diagnoses and their associated candidate models, the 
next phase of the diagnosis process is quantitative refinement of the qualitative can- 
didate diagnoses and their associated models through parameter estimation and data 
fitting, followed by tracking of the fit of subsequent observations to the candidate mod- 
els. The goal is to at least provide a probabilistic ranking of foe plausible candidates, if 
not a unique model (and hence diagnosis). 


2 We may cut off back-propagation along foe mode sequence beyond a tune limit. 


As observed in the previous section, die model associated with the candidate qualita- 
tive diagnosis, Mode is underconstrained. Both die time of the fault mode occurrence, 
/ r and die parameters associated with the faulty behavior Op are represented as ranges 
and must be estimated. Further, the candidate qualitative diagnoses were generated from 
initial observations of aberrant behavior, and their consistency can be further evaluated 
by monitoring the qualitative transients associated with each candidate. The refinement 
process is performed by a set of trackers [21], one for each candidate diagnosis and 
associated model. Each tracker comprises both a qualitative transient analysis compo- 
nent and a quantitative model estimation , component. The two components operate in 
parallel as described below. 

Qualitative Transient Analysis 

The qualitative transient analysis component performs a further qualitative analysis of 
the consistency of candidate qualitative diagnoses based on monitoring of higher-order 
transients whose manifestation is seen over a longer period of time. If the transients 
of a candidate qualitative diagnosis do not remain consistent with subsequent observa- 
tions, fae candidate diagnosis will be eliminated and the model estimation component 
informed. The technique we employ is derived from techniques for qualitative monitor- 
ing of continuous systems. Details of the algorithm appear in [18]. 

Model Estimation 

The purpose of the model estimation component is to perform quantitative model fit- 
ting, i.e., to provide a quantitative estimate of the parameters of the models and to assign 
a probability to each of the candidate models (and hence candidate diagnoses), given 
the noisy observed data. In particular, given a candidate model, M ode the model es- 
timation component uses parameter estimation techniques to estimate both the time at 
which the failure occurred, l f, and the value for the parameters, Op, associated with the 
conjectured failure mode. In this paper we discuss two alternate approaches to our time 
and parameter estimation problem. The first approach is based on Expectation Maxi- 
mization (EM) (e.g., [8]), an iterative technique that converges to an optimal value for 
If and Op simultaneously. The second approach we consider employs General Likeli- 
hood Ratio (GLR) techniques (e.g., [5]) to estimate the time of failure / p, and then uses 
the observations obtained after the failure to estimate the fault parameters, 6 p, by a least 
squares method. As described in Section 3, the outcome of both approaches is a unique 
value for Ip and Op and a measure of die likelihood of Mode given the observations. 
The proposed approaches to model fitting have trade-offs and we are currently assess- 
ing the efficacy of these and other alternative approaches through experimentation. 

EM-Based Approach The Expectation Maximization (EM) algorithm (e.g., [8]) pro- 
vides a technique for finding the maximum-likelihood estimate of the parameters of an 
underlying distribution from a given set of data, when that data is incomplete or has 
missing values. The parameter estimation problem we address in this paper is a vari- 
ant of the motion segmentation problem described in [24]. Here, we define die basic 
algorithm and the intuition behind our approach. (See [8] for more details.) 

The time of failure, I p = [//, / „] of our candidate qualitative diagnosisdictates the 
mode in which the failure is conjectured to have occurred. Let us call this mode // *. 
The behavior of our hybrid system in mode /<, is described by die continuous function 


f /li9 with known parameters At some (to be estimated) time point ij? within the 
predicted time period of /*,, we have conjectured that die system experienced a fault 
which transitions it into mode ftp. The behavior of our hybrid system in mode ftp is 
described by the continuous function / /tP , with unknown parameters, Op. We also have 
a set of data points O = [x^ (//),... ,**/,*(/«)] Q O, which either reflect the behavior 
of the system under f /ti or under f ftr . 

Given all this information, our task is to find l) values for parameters Op , and 2) an 
assignment of the data points O to either }tj or ftp so that we maximize the fit of die 
data to the two functions. The assignment of data points will in turn tell us the value 
of Ip. EM provides an iterative algorithm which converges to provide a maximum- 
likelihood estimate for Op given O , i.e., roughly we are calculating the likelihood of 6, 
L(0) = P(0 | Op, Mode). 

The basic EM algorithm comprises two steps: an Expectation Step (E Step), and a 
Maximization Step (M Step) [24]: 

• Select an initial (random) value for Op. 

• Iterate until convergence: 

- E Step: assign data points to either f fli (9i ) or which ever fits it best. 

- M Step: re-estimate Op using the data points assigned to f fl P (Op). 

The assignment of data points to //; and ftp provides an estimate for Ip, We may 
exploit the fact that the assignment of data points is temporally correlated with all points 
before I/ belonging to //<, and all points after If belonging to ft /. We may also exploit 
the fact that data points at the beginning of the interval will belong to /< ;, while those 
at the end will belong to ftp. These task-specific qualities help our algorithm converge 
more quickly. 

EM provides a rich algorithm for maximum-likelihood parameter estimation when 
we don’t know the value of/p. In some hybrid diagnosis applications, depending upon 
the sensors in our system, and the level of noise in die sensors, we may be able to de- 
velop monitoring techniques that will help isolate a reasonable value for Ip, minimizing 
the need for iteration in EM. In such cases, an alternative to the EM-based approach is 
to first estimate Ip using the Generalized Likelihood Ratio (GLR) method [5], followed 
by parameter estimation of Op. 

GLR + Least Squares Approach Here, we divide the parameter estimation problem 
into two parts: (i) estimate the time of failure, Ip, using die Generalized Likelihood 
Ratio (GLR) method, and (h) apply a standard least squares method for parameter esti- 
mation. The intuition is that solving the problem in two parts simplifies the estimation 
process, and very likely mitigates die numerical convergence problems that arise in 
dealing with complex higher-order models. 

The GLR method for detecting abrupt changes in continuous signals is described 
in [5]. We have applied it to fault transients analysis in complex fluid thermal systems 
[16]. Here we provide an overview of the method for die single parameter case. Assume 
that the signal under scrutiny is a time-indexed sequence of random variables y(k), with 
probability density function, po, (y) in desired mode /./,*, and po r (y) in fault mode ftp. 
y is either contained in * 0 hs or computed from We assume that a fault causes an 
abrupt change in p(Jfc). In the case of the AERCam, y captures the difference between 
the observed and expected values of the, e.g., acceleration, as predicted by the model. 


The central quantity in the change detection algorithm is the cumulative sum of the 
log-likelihood ratio for a window of observations between times in and n, 
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Again, this ratio is a function of two unknowns: Ip and Op. The common statistical 
solution is to use maximum likelihood estimates for these two parameters, resulting in 
a double maximization: 


fin = max sup S;;,(6»r). 
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If we assume that probability density functions, po ; (y) and po r (y) are Gaussian, 
then g n reduces to: 
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where u and <r‘f are the mean and variance for p 0/ (y), respectively. 

When processing a sequence of samples, the point of abrupt change, / p, is computed 
from mln{n : y n > //.}, where h is an appropriately defined threshold. Hence, the 
smaller the value of h, the more sensitive the function to change, and unfortunately to 
false alarms, so h must be set carefully. 

Once / p is estimated, data points observed after / p, are used to estimate the parame- 
ter, Of? for a hypothesized fault using regression techniques. In the case of the AERCam, 
the position vector of the AERCam is modeled as a set of quadratic functions in terms 
of the thruster force. These functions contain one unknown, 0 p, the parameter that cor- 
responds to the degree of degradation in the faulty thruster. The least squares estimate 
for Op is computed, and the the measure of fit of the candidate model to the observed 
data used to estimated the probability of the candidate model (and hence, diagnosis). 

Model Comparison 

From the model estimation component, each tracker computes the likelihood of its 
model Mode, and hence of the associated candidate diagnosis {C.//p. / p. Op), as a 
measure of fit of the observations to the model. As new data are observed. Op 

and /p, are adjusted and P{Modc | computed. If the likelihood of Mode 

falls below a predefined acceptable likelihood threshold, a, then its tracker is termi- 
nated, and the associated candidate diagnosis (C ? /tp ? Ip. Op) removed from the list of 
candidate diagnoses. Tracking terminates when a unique diagnosis is obtained, or when 
the diagnoses are sufficiently discriminated to determine suitable controller actions. 

5 Related Work 

The specific problem of diagnosing hybrid systems has received little attention to date, 
although there is much related woik. Within the AI community, there has been a great 
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deal of research on diagnosing static systems (e.g., [14]), while much less on diag- 
nosing discrete dynamical systems (e.g., [17,25]), and qualitative representations of 
continuous systems (e.g.,, [18]). Within the FDI community, the largest proportion of 
research has focused on diagnosing continuous systems (e.g., [13, 1 1]). The most com- 
mon model-based approaches use observer schemes(e.g., [ 12, 20]), where the goal is to 
design residual generators based on observed discrepancies, such that individual resid- 
uals are sensitive to a particular subset of faults. There is also complementary work by 
Basseville [4], using model-based statistical processing techniques for early fault de- 
tection and residual identification. [18] perform residual generation and analysis task in 
a qualitative framework to address some of die computational issues that arise in han- 
dling the complex dynamics that occur in fault transients, with some preliminary work 
on building multiple observers for hybrid systems [19]. Diagnosis of discrete-event sys- 
tems has also been studied within the FDI community (e.g, [22, 15]). Fabre et al. [ 10] 
have employed stochastic Petri nets based on a Hidden Markov Model probabilistic 
scheme for alarm analysis. Unfortunately, it is not clear how to systematically derive 
such representations from die physical system models that we work with. 


6 Summary 

In this paper we addressed the problem of diagnosing hybrid systems. The main con- 
tributions of the paper are 1) formulation of die hybrid diagnosis problem as model 
selection; 2) the exploitation of techniques for qualitative diagnosis of continuous sys- 
tems to reduce the diagnosis search space; and 3) the use of parameter estimation and 
data fitting techniques for evaluation and comparison of candidate diagnoses. This work 
continues with experimental analysis of the proposed techniques, and a more formal 
characterization of our approach in terms of Bayesian model selection. 
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